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ABSTRACT 


THE EFFECT OF ELASTIC BOUNDARY CONDITIONS ON THE 
DYNAMIC RESPONSE OF RECTANGULAR PLATES 

Terry K. Brewer 
Old Dominion University 
Director: Dr. Chuh Mei 

Natural frequencies and forced steady-state harmonic 
response for the vibration of uniform rectangular plates with 
edges elastically restrained against rotation and transverse 
translation are addressed. A single mode Rayleigh-Ritz 
solution is derived using functions that describe the normal 
inodes of vibration of a beam whose ends are elastically 
restrained. The finite element solution is obtained for 
comparison. MACSYMA symbolic manipulation system is imple- 
mented as an aid to the mathematical rigor of the Ritz 
approach and NASTRAN finite element code is used to model 
the mechanical system. Comparisons are made to published 
results and the solutions of this study are found to give 
lower frequencies for some values of boundary restraint. 
Steady-state harmonic amplitudes of displacement and acceler- 
ation are found to agree favorably for the two solutions. 

Low predictions of steady-state strain from NASTRAN result 
in some cases when compared to the Ritz values. Finally, a 
subjective assessment is made about the merit of using 
MACSYMA and NASTRAN. 
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CHAPTER 1 
Introduction 

Research studies of plate vibrations under various 
loadings, plate geometries and panel types have received 
wide attention since the 1950's due mainly to the appli- 
cations of these studies to the fatigue failures of air- 
craft structures. Until recently, design of the structures 
has been largely based on empirical nomographs (1-3) but as 
new materials (including composites and high temperature 
materials) gain wider applications to aercspace vehicles, 
the development of design charts is economically limited. 
Currently, these new materials are generally reserved for 
"low risk" applications. As a result, there is a need to 
produce and improve analytical procedures to accurately 
predict panel response; 

Since the predicted and measured bending strain of 
conventional aluminum rectangular panels can typically 
differ by a factor of two or more, studies have been made 
to assess the accuracy of classical linear small deflection 
plate theory. Roussos, Heitman and Rucker (4) encountered an 
even greater discrepancy in 1986 as predicted strain values 
were approximately three times higher than those measured. 
Although these results may be "safe" from a design standpoint. 
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the optimization of design is sacrificed. This large bias 
error prompted the investigation of the effect of boundary 
conditions on the response of rectangular plates (5). In part- 
icular, the interest was focused on edges that are restrained 
by both rotational and translational boundary springs. 

Previous experimentation was conducted on four panel 
types including isotropic (aluminum), orthotropic and two 
anisotropic plates. The mounting of the boundaries was done 
utilizing a rubber gasket scheme on both sides of the panels 
with double gasket thickness on the load source side. The 
load source delivered an acoustic excitation at levels of 
80 to 120 dB. These low levels of normally incident plane 
waves were used to apply a uniform spatial load at constant 
temperature. The horn speaker produced low frequency (60 to 
2000 Hz.) noise as the sound field was mapped and the pressure 
spectrum was analyzed. The panels were measured for fundamental 
resonance frequencies and critical damping ratio. However, 
the edge restraint values were not experimentally determined 
but were considered uniform on each edge by virtue of the 
rubber gaskets. 

The analytical models employed for comparison of measured 
acceleration and strain response displayed only rotational 
boundary restraints with edges rigidly supported in the tran- 
sverse direction. The accompanying analysis was performed by a 
Rayleigh-Ritz frequency response formulation and a supplemental 
response calculation was done on the isotropic panel -by a 
finite element method. The boundary values used were 
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determined by requiring good fundamental frequency prediction 
and close agreement to ratios of x and y strain measurements. 
The measured accelerations closely matched those predicted. 
Unfortunately, the measured and predicted strain response 
differed by the aforementioned bias for all panels tested. 

In a search for causes, the finite element model was extended 
to include the elastic transverse boundary conditions. In the 
comparison to measurements, the relaxation of transverse edge 
motion accompanying rotational edge elasticity could not 
account for such a large bias. Fault for the discrepancy was 
tentatively placed on the strain gages chosen for testing. 

To date, most of tne analysis of elastically restrained 
plates has been done using the Rayleigh-Ritz energy form- 
ulation although the Galerkin method and numerical techniques 
have been utilized (6,7) . In the Rayleigh-Ritz approach, the 
maximum strain energy is equated to the maximum kinetic 
energy of the plate and by assuming a finite linear series of 
displacement functions that satisfy the geometric boundary 
conditions , the natural frequencies, mode shapes and 
consequently, the overall response may then be solved for. 

Much of this work has been presented using polynomials as the 
assumed displacement functions since the resulting calculations 
and computations are made simpler. This approach has been 
shown to be adequate to obtain natural frequencies and mode 
shapes for lower modes of vibration (8). It is also advan- 
tageous when the study includes nonrectangular plates, non- 
uniform thickness and plates with geometric discontinuity 
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because of the simplicity of the assumed polynomial functional 
relation (9,10). 

Since the accuracy of the Rayleign-Ritz approximation 
depends, to a large extent, on the "closeness" of the assumed 
displacement functions to the true mode shape, many analysts 
feel the use of characteristic modes of beam vibration in the 
method is a better choice. The underlying assumption is that 
in a given coordinate direction, a slim strip of the plate 
will approximate the behavior of a vibrating beam. The plate 
is then modeled as a continuous collection of beams in two 
directions. This choice of functions has been demonstrated 
to give greater accuracy especially for the natural frequencies 
of higher modes and in the case of mixed classical boundary 
conditions (11). Classical boundary conditions are traditional- 
ly regarded as clamped, simply supported and free edges and can 
be mixed among the four edges to yield 21 combinations. 

It should be noted that the Rayleign-Ritz method and 
other techniques used to solve problems of rectangular plate 
vibrations are merely approximations, in general. Exact 
solutions to the differential equation of transverse 
vibrations exist only for the cases of paired simply supported 
parallel edges. Of the 21 possible combinations of classical 
boundary conditions, these exact solutions comprise only six 
cases. By excluding the free edge cases and imposing mounting 
constraints, the number of classical boundary exact solutions 
drops to three. These three are the cases of a plate with all 
boundaries simply supported, three edges simply supported with 
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one end clamped and two opposite edges simply supported with 
clamped conditions on the remaining two boundaries. Since, in 
reality, there are no perfectly clamped or simply supported 
edges, it is reasonable to restrain the boundaries with 
continuous rotational and/or translational springs. The 
classical edge conditions become special cases with approp- 
riately assigned spring constants. For example, the simply 
supported edge is achieved by imposing a rotational spring 
value of zero and a transverse foundation spring value of 
infinity. With these considerations included, practical 
mounting of a rectangular plate is more realistically modeled. 

In 1983, Warburton and Edney (12) presented the Rayleigh- 
Ritz method utilizing suitable combinations of modal functions 
of beams with classical boundary conditions to study elastically 
restrained plates. As an example, the natural frequencies of a 
panel with parallel edges restrained against rotation are inter- 
mediate between those of plate with opposite edges simply 
supported and those with opposing edges clamped. Tne assumed 
mode shape is a summation of simply supported-simply supported, 
simply supported-c lamped , clamped-simply supported and clamped- 
clamped characteristic beam functions in the appropriate 
coordinate direction. In this way, combinations of the limiting 
values of the rotational restraint parameters on each edge 
are included in the analysis and the classical functions can 
be weighted to give the intermediate values. By incorporating 
the free edge and sliding edge conditions of beams into the 
procedure, translational boundary restraints may be included 
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also. This method proved to give good accuracy when compared to 
the results of previous authors and was shown to be relatively 
easy to use since the integrals of beam functions that result 
in the Rayleigh-Ritz approach and the characteristic frequency 
parameters for classical boundary conditions are well docu- 
mented. The underlying purpose of Warburton and Edney was to 
give designers rapid and inexpensive estimates of natural 
frequencies of practical plates with applications to many modes 
of vibrations, aspect ratios and edge restraint parameters. 

One of the methods used for comparison by Warburton and 
Edney was given by Carmichael in 1959 (13). Carmichael also 
used the Rayleigh-Ritz procedure. The analysis treated plates 
with rotational edge restraint only whereby all boundaries were 
considered rigidly restrained against translation. Accordingly, 
Carmichael assumed displacement functions that represent the 
normal modes of vibration of a beam whose ends are elastically 
restrained against rotation and rigid in the transverse direc- 
tion. This approach necessitated the development of integrals 
of these elastically supported beam functions and he also 
generated a closed form approximate frequency expression from 
a single term solution. This approximate frequency was shown 
to be quite accurate as compared to the frequency from a 
multi-term solution. 

Other studies of elastically restrained plate vibrations 
are usually compared to those mentioned when relative merit 
is sought. Analysts have given quick methods for natural freq- 
uency calculations and tables of values for the use of a 
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given approach. However, there is a loss of generality assoc- 
iated with some derivations as the authors restrict themselves 
to a specific problem with simplified assumptions. For 
example, only recently has attention been given to the case 
of edges possessing two distinct types of boundary springs. 
Prior to this, a presentation might treat problems exhibiting 
only symmetric rotational boundary springs and/or square 
plate geometries. In order to further these methods to a 
larger class of problems, it becomes necessary to "rework" 
the entire problem. Thus, the application of an approach 
can be long and arduous. 

Another difficulty associated with the implementation 
of an analytic technique is the brevity of published treatise. 
Often, the answers sought and presented are limited to the 
extraction of natural frequencies and mode shapes related to 
free vibration of the plate so that examples of numerical 
results of displacements, accelerations and strains are not 
implied or referred to. As a result, the engineer using the 
technique to predict forced response must finish the problem. 
This is not always simple if functional expressions for 
numerical calculations are desired. The mathematical manipu- 
lation required to obtain a completed solution can be algeb- 
raically difficult and seemingly impossible. Usually, an 
assessment is made to decide if the extra analysis is practical 
in both time and expense of the work. If the labor is pursued, 
it becomes necessary to evaluate the results by comparison to 
real problems. 



Brewer 


8 


The comparisons made by many authors is a third pertinent 
topic. Quite often, the evaluation of an operational algorithm 
is limited to a comparison with a procedure that is similar in 
derivation. Intuitively, the two approaches yield similar 
results. The studies are not usually compared to dissimilar 
formulations of the same problem. Since prediction of 
structural behavior is the ultimate goal of these studies, 
comparisons are required to different solutions and to 
measured response from experimentation. 

In order to adequately address the topics of elastically 
restrained rectangular plate vibrations, the present treatise 
again utilizes the Ray leigh-Ritz method. However, the 
assumed displacement functions are given by tne normal modes 
of beams whose ends are unsymmetrically restrained against 
rotation and translation. This is an extension of Carmichael's 
approach and necessitates the derivation of the beam solution 
and the generation of the resulting beam integrals. In this 
way, unlimited combinations of distinct rotational and linear 
spring parameters can be included. Variation of aspect 
ratio can also be considered. 

Because of the complexities in the mathematical 
rigor associated with this type of solution, a symbolic 
manipulation program is used as an aid in the analytic 
derivation. The program chosen is MACSYMA since it is very 
comprehensive in it's mathematical capabilities and it is 
well documented. The differentiation and integration required 
are effectively executed as program commands. The evaluation 
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of boundary condition:; and simultaneous equation algorithms 
are available to the user. It is possible to generate an 
analytic solution and to program numerical procedures of 
problems in rectangular plate vibrations completely on the 
computer. This has great potential when compared to "hand- 
cranking" an analytical derivation. 

The comparisons of solutions includes results from 
previous authors where applicable as well as the comparison 
to a finite element formulation. The well Known finite 
element code, NASTRAN,, is utilized. The choice of NASTR^N 
was also based on availability and comprehensive capabilities. 
These comparisons allow the assessments of accuracy of 
solutions and numerical results. This evaluation also gives 
an estimation to the relative contribution value of 
available engineering and mathematical computer software. 

More specif ically, the purpose of this thesis is to study 
the effect of elastic boundary conditions on the dynamic 
displacement and strain frequency response of rectangular 
uniform plates. The loading is assumed to have a deterministic 
frequency spectrum that can be represented by a finite 
series of harmonic inputs. The loading amplitude is assumed 
uniform over the plate surface. The direction is to normal 
incidence. Parameter studies are made for plates at several 
aspect ratios for predicted natural frequencies versus elastic 
restraint values. Studies include trend analyses for steady- 
state displacement and strain response. 



Brewer 10 


Two methods are used to conduct these studies. They are 
the Raylei (jh-Ritz method and the finite element method. 
Comparisons are made to previous accepted results to verify 
the two methods as solution techniques. The Rayleigh-Ritz 
approach is aided by MACSYMA symbolic manipulation program. 
NASTRAN finite element code is utilized. A subjective assess- 
ment is made of the relative contribution of MACSYMA and 
NASTRAN as computational aids. 

The contributions of this thesis include the use of 
beam functions representative of the elastic boundary con- 
ditions in the Rayleigh-Ritz approach and the extension of 
the method to predict measurable response. 

Classical linear small deflection plate theory is 
assumed throughout the analysis. 
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CHAPTER 2 

Rayleigh-Ritz Energy Formulation 
Theoretica 1 Basis 

In 1877, Lord Rayleigh used the principle of conservation 
of energy to obtain an expression for the fundamental natural 
frequency of an undamped linear structure under free vib- 
ration (14). He argued that as the structure vibrates har- 
monically at it's natural frequency, the potential energy 
given when vibration is at maximum amplitudes is exchanged for 
the energy of a purely kinetic state when amplitudes are zero. 
The maximum potential energy of bending, V max b , of a plate 
is determined by the strain energy stored in the plate. It 
should be noted, however, that V max , b from plate bending is 
not necessarily equal to the total maximum potential energy, 
V max , since additional terms in V max may result from other 
energy storage sources such as boundary springs, for example. 
For the case of the aforementioned classical boundary con- 
ditions, v max>b = v max . 

The maximum kinetic energy, T max , is taken while the 
plate vibrates at it's natural frequency and the fundamental 
frequency can be calculated by minimizing the following 
equation , 

T max " v max = stationary value (1) 
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Ideally, the expression given in equation (1) is identically 
zero. However, this only occurs when tne true mode displace- 
ment function is used to determine V max and T max . Since 
exact solutions to rectangular plate vibrations are not 
known, in general, equation (1) is minimized with respect 
to the assumed displacement functions. This reduces the 
resulting nonzero error to a minimal stationary value. If 
the assumed displacement function satisfies all the boundary 
conditions, the calculated natural frequency obtained has 
been shown to be an upper bound. That is, the true structural 
frequency is lower than that calculated. The reason for this 
is that using an assumed mode in the procedure is equivalent 
to an additional constraint which raises the potential 
energy and calculated frequency. The Rayleigh formulation 
frequency approaches the exact from above as the greater ac- 
curacy criteria dictates. 

For example, the bending strain of a uniform rectangular 

* 

plate can be shown to be given by, 


V 

max 



+2 M W 'xx w 'yy 2 

+2 ( 1- yu. ) W, xy 


)dydx 


(2) 


where for now it is assumed V max>b = V max , and 

where the comma denotes partial differentiation with respect 

to the variables following the comma. 


* Notation convention is given in Appendix A. 
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The maximum kinetic energy is given by, 


T 


max 



W (x,y) dydx 


(3) 


where w = w(x,y)e is the assumed transverse displacement 

function over the spatial domain vibrating at the fundamental 
frequency, U) . The i denotes^ -1 ) . 

For the ideal case (stationary val!ue=0), substituting 
equations (2) and (3) into (1) and solving yields, 

V„ 


2 

CO = 


max 



W(x,y) dydx 


(4) 


Equation (4) is called the Rayleigh Quotient. 

In order to obtain frequencies beyond the fundamental, 

Ritz developed an extension of the Rayleigh formulation. The 
estimated mode function is assumed to be a finite linear series 
with arbitrary constant coefficients, C r . The minimization of 
equation (1) then becomes, 

( V - T ) = 0 , ( r=l,2,3,...n) 

(5) 


a 


Equation (5) results in n homogeneous equations that can be 
solved for n-1 constants in terms of the remaining constant. 
When the finite series is truncated at r=l , equation (5) 
reduces exactly to equation (4) for the determination of 
fundamental frequency while C} is arbitrary. The remaining 
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constant represents an amplitude and for free vibration of 
the plate, it is dependent on the initial conditions. For the 
case of steady state harmonic forced response, the amplitude 
constant is dependent on the amplitude and frequency of the 
forcing function. When n>l, equation (5) can be written in 
matrix form and effectively becomes an eigenvalue problem 
with mass and stiffness contributions arising from the 
kinetic and potential energies, respectively. The natural 
frequencies are obtained from the condition that the deter- 
minant of the system of equations is zero. 

C onsiderations of Boundar y Conditions 

Ideally, all boundary conditions are satisfied by the 
chosen displacement functions. However, previous methodology 
(14) has shown that sufficient accuracy of frequency prediction 
can be achieved by satisfying only the geometric boundary 
conditions (ie. displacement and slope) while neglecting 
the dynamic or natural boundary conditions ( ie . shear and 
moment). The accuracy is dependent on the true type of edge 
conditions for a case under study. For example, the clamped 
edge has only geometric boundary conditions, W(o,y)=0 and 
w, x (o,y)=0 , and the Rayleigh-Ritz method gives quite good 
frequency results with appropriately chosen functions. When 
the edges are restrained elastically, the boundary conditions 
can be considered natural only. The choice of functions that 
satisfy geometric conditions can result in loss of accuracy. 

For the plate shown in figure (1), the rotational 
boundary conditions are given by, 
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Figure (1). Mechanical System Under Study 
R 1 W'x* 0 '^ D< W 'xx (o 'y> + H- W »yy(°»y)) 


R 2 W, x (a ,y ) == "D(W, xx (a, y ) + /X W, (a,y)) 


R 3 W, y (x,o) = D(W, yy (x,o) + W » xx ( x , o ) ) 


R 4 W, (x,b) -D(W, yy (x,b) + /i w, xx (x,b) ) 


(6A) 


(6B) 


(60 


(6D) 


where R^ is the bending moment (per unit length) for unit 
rotation on the i th edge, (i=l,2,3,4). 

The translational boundary conditions of the plate are 


given by. 
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Kj W(o,y) = - D(W, xxx (o,y) + (2- yx )W, xyy (o,y)) 
K 2 W(a,y) = D(W, xxx (a,y) + (2- /X )W, xyy (a, y ) ) 
K 3 W(x,o) = -D(W, yyy (x,o) + (2- yx )W, yxx (x,o) ) 


(7A) 

(7B) 


(70 


K 4 W(x,b) = D( W, yyy ( x , b ) + (2- yx )W, yxx (x,b) ) 


( 7D) 


where is the transverse edge reaction (per unit length) 

for unit displacement on the i th edge, (i=l,2,3,4). 

Equations (6)-(7) are simplified if in addition to 
uniform elastic stiffness, twisting on each edge is not 
restrained. That is, the moments about an axis parallel to 
a given edge and the vertical shear on that edge are directly 
restrained by the boundary springs. The boundary conditions 
become , 

R 1 ^'x' 


l 2 w 'x 


‘3 w ' y 


<4 W/y 


(o,y) 

= D ( w , xx ( ° » y ) ) 

( 6A ' ) 

(a,y ) 

= -D( W, xx (a , y ) ) 

(6B') 

(x,o) 

= D(W,yy (x,o) ) 

( 6C ' ) 

(x, b) 

= -D(W, yy (X,b) ) 

( 6D ' ) 

*y) = 

_D(w 'xxx (o 'y ) ] 

( 7A ” ) 

»y> = 

D( w , xxx (a,y)) 



(7B') 
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Kj W ( X , O ) = -D(W, yyy(X r O) ) 


( 7C ' ) 


K 4 W ( X , iD ) = D(W, yyy (X,b) ) 


( 7D ' ) 


The boundary conditions given by equations (6')-(7') 
are assumed in this analysis and the beam functions used 
are chosen so as to satisfy these conditions. All in-plane 
motion at the boundaries is prohibited since the assumption 
of linear small deflection theory renders these motions 
negligible. That is, for small amplitude plate vibrations, the 
transverse bending and the in-plane membrane oscillations 
are uncoupled and independent. With this in mind, transl- 
ational motion at the edges is meant as transverse motion 
only . 

Choice of Assumed Displacement Functions 

The general form of the finite series to represent 
displacement in the Rayleigh-Ritz method for harmonic steady 
state rectangular plate vibrations is given by, 


i CJ t 

w(x,y ,t ) = W(x,y ) e 


( 8 ) 


where w(x,y,t) is assumed to be separable in x,y and t. The 
spatial variation, W(x,y) , is, 



m=l n=l 


^mn X m^ x ^ Y n^^ 


( 9 ) 


such that X m (x) and Y n ( y ) are chosen to satisfy the 
boundary conditions while m and n are modal identification 
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indices. The beam functions selected are given by, 

CL n\ x Cl 

X m (x) = B m sinh ( — ci) + C m sin ( — a ) + 


a m x Q m x 

D m cosh( — a — ) + cost — a — ) 


Y n (y) = E n sinht 


.0 n y P ny 

d ) + F n sinCTi ) + 

R n y £ n y 

^"d ) + cost : 


G n cosht ' d ) + cosC“E ) 


( 10 ) 


( 11 ) 


where B,C,D ,Cl,E,F,G and p are constants 
obtained from the application of equations (6')-(7'). The 
results are given in Appendix B. 

By noting the similarity of x to y in the equations (10)- 
(11) and also in equations (6')-(7'), the respective constants 
are of like form. That is, once B , C , D and d are deter- 
mined from equations (6A')» (6B'), (7A'),(7B') and (10), E , 

F , G and p result by substituting y for x, E for B, F for 
C, G for D, p for cl for m, Rj for , R 4 for 
R 2 , Kj for Kj , K 4 for K 2 and b for a. 

D erivation of^ Solution 

The differential equation that governs the motion of a 
uniform plate subject to an external loading is given by, 



p h w + c w = p 


( 12 ) 


where p h is mass per unit area ( p= mass density) and c is 
the damping coefficient. 

The pressure, p = p(x,y,t) , is assumed to be harmonic 
and uniformly distributed which results in an harmonic 
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displacement, w = w(x,y,t). Using complex notation, 

i (jJ t i U) t 

p = p(x,y,t) = P(x,y) e = P e 

and, i OJ t 

w = w(x,y,t) * W(x,y) e 


(13) 

(14) 


where P(x,y) and W(x,y) are the spatial amplitudes of 
pressure and displacement, respectively. For a uniform 
pressure, P(x,y) = P . 

Since none of the Boundaries are simply supported, the 
approximate solution comes from an application of the Rayleigh- 
Ritz approach. This method in dynamic response problems 
involves forming the kinetic energy, T, the potential energy, 

V, and the work done,U, by the pressure on the plate. The 
solution for W(x,y) is pursued by minimizing the following 
equation , 

T - (V + Q)= stationary value (15) 

The kinetic energy is calculated by equation (3). 

The potential energy is formed by summing the contri- 
butions from the plate bending strain energy and the strain 
energy in the boundary springs. The bending strain energy 
can be determined from equation (2). The potential energies 
associated with the rotational elastic boundary conditions 
are. 
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(16) 


The strain energies associated with the transverse boundary 
springs are, 



where all spring coefficients ( 
assumed constant along respective edges. 
The total potential energy can then 

V - v R + v K + v maXyb 

The work done by a uniform pressure 


2 

W (a,y ) dy 


* 2 

W (x,b) dx 

(17) 

( i=l, 2,3,4)) are 
be written as, 

(18) 

field, P , is given 


as , 



dy dx 


(19) 


The choice of W(x,y) was given in equation (9) so that 
the minimization of equation (15) becomes, 



8v_= 

8 Vi 


<3q 

3 Anm 


( 20 ) 


Since the response of the plate is ruled by the natural 
frequencies and geometric mode shapes, the calculation of 
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natural frequencies is required.. This is done by forming the 
Rayleigh Quotient of equation (4) so that the radian frequency 
of the mn th mode is given by. 



( 21 ) 


where equation (21) can be written as a matrix which is a 
system of simultaneous equations. The natural frequencies are 
obtained when the determinant is set to zero. It involves the 
evaluation of many difficult integrals. For this reason, 

MACSYMA is very helpful. The evaluation is given in Appendix C. 

Combining equations (20) and (21), the minimization 
problem becomes, 

_dr_( 1 - ^ron ) ~ 3 q 

3 A,n W 2 d (22) 


whore equation (22) is also a matrix of simultaneous equations 
which allows the calculation of the constants, A^, by 
substitution of equations (3), (9)— ( 11) and (19). These 
results are described in Appendix D. 

The calculation of the coefficients, h mn , at this point 
makes the solution for the spatial distribution of displacement, 
W(x,y) , complete. 

The resulting strains are then calculated by the small 
deflection theory as follows. 
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and , 


i2 TT ft 

€ x = 2 w 'xx = “ 2 W 'XX e 


i2 TT ft 

€ y = - z w , yy = - z W, yy e 


(23) 


(24) 


where 2 TT f = CJ . The maximum strains are calculated at the 
top and bottom surfaces of the panel at z=+h/2. 

The magnitude of mean-square strains at a frequency, f, 
are determined by. 


and , 



= -U € € * 

2 x x 




= -I- € €* 

2 y y 



(25) 


( 26 ) 


where * denotes the complex conjugate. 
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CHAPTER 3 

Symbolic Manipulation as an Aid to Solution 

* 

Introduction and Description of MACSYMA 

Historically, machine aided mathematics has experienced 
an arithmetical development. The focus has been to provide a 
means to process large amounts of information numerically. 
This is reflected in the widely held desire to statistically 
evaluate numerical data. As a result, the majority of soft- 
ware systems available are designed for digital data proces- 
sing. Any analytical relation related to the numbers is 
usually derived prior to computer use and incorporated into 
the data processing or the relation is discovered with fur- 
ther numerical evaluation such as curve fitting. The computer 
has proven to be very capable to perform statistical analysis 
but the accuracy and efficiency are not always good. Large 
sets of data require large storage space in the machine and 
truncation with roundoff is not uncommon. 

Since the computer combines and arranges numerical 
quantities based on mcithematical expressions that manipulate 
variables, parameters or any general symbol, it is not 
unreasonable to expect, these same mathematical expressions 


* Acronym for MAC 's symbolic Manipulation System 
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to bo derivable on the machine. During the execution of an 
algorithm, the computer need not discriminate between num- 
bers and symbols until the user so specifies with an output 
request. Research and development of computer aided mathematic 
is best exemplified by MACSYMA. 

MACSYMA 's initial development was done under the 
direction of Professor Joel Moses at Massachusetts institute 
of Technology in the latter part of the 1960 's (15). Since 
then, it has become commercially available through Symbolics, 
Inc. The history of MACSYMA 's development shows continuous 
contributions to improvements and system maintenance by many 
authors. The system has evolved into a large and sophisticated 
symbolic manipulation software package that gives the user 
access to many otherwise unavailable mathematical techniques. 
This has produced a library of user shared algorithms that 
have been employed in such fields as experimental mathematics, 
acoustics, celestial mechanics, computer aided design, 
structural mechanics and numerical analysis. The types of 
problems that have been addressed include spectral analysis, 
helicopter blade motion, finite element analysis and plate 
vibrations, among many others (16). 

MACSYMA is a lisp program (recently estimated at 
300,000 lines of lisp code (16)) that has been catalogued 
to include basic algorithms and axioms of matnematical tneory. 
The growth of MACSYMA is not unlike the evolution of math 
itself in that advanced theory is dependent on the validity 
of fundamental assertions. This " building block " process 
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is apparent in mathematical proofs or derivations whereby 
any conclusion must logically follow preceding steps of the 
argument and MACSYMA has experienced this same " building 
block development. The program has attained a level of sop- 
histication that offers on-screen two-dimensional represent- 
ations (very similar to handwritten form) of advanced cal- 
culus and linear algebra. Accordingly, the program commands 
for the manipulation of user defined variable expressions 
are similar to the mathematical functions desired. For 
example, in order to factor an expression, the MACSYMA 
syntax is, 

factor (expression) ; 

or to expand a rational expression by cancelling a common 
divisor, multiplying out products of sums, etc, the syntax 
is , 

ratexpand(expression) ; 

or to simplify a trigonometric expression by the implemen- 
tation of the identity, sin 2 (x) + cos 2 (x) = 1, the command 
is , 

trigsimp( expression ) ; 

Also, the arithmetical operators used to enter and create 
equations are the same as those used in many program languages 
where ** = exponentiation, / = division, etc. As a 
result, it is not difficult to begin the application of 
MACSYMA capabilities and subsequent requests are user dis- 
covered in a friendly manner. 
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Perhaps the most valuable asset given by MACSYMA and 
other available symbolic manipulation software is the release 
from computational details of problem solutions. Many solution 
derivations involve long and difficult intermediate calcu- 
lations where a strong chance for error is present if the 
computations are done by hand. Even commonly used tables of 
integrals contain errors. A computer algebra system puts the 
mathematics of solution in the hands of the practicioner . 

The answers returned by MACSYMA are given in general form 
according to the degree of generality requested. This frees 
an engineer to directly define and apply general solution 
techniques and algorithms to a class of problems and allows 
the examination and development of theory from the definitions. 
The insight to be gained from engineering theory and physical 
phenomenon is not lost in tne details of derivation. Since 
mathematics is a tool of the engineer, it is advantageous to 
let the computer assume the role of the analytic work-horse. 
Examples of MACSYMA Capabilities 

The capabilities of MACSYMA are too numerous to mention 
all here. Examples, with an emphasis on those related to the 
current study, include the ability to perform definite and 
indefinite integration, differentiation, evaluation of limits 
and series summations, trigonometric and hyperoolic function 
manipulation as well as matrix and tensor manipulation. The 
program can also evaluate expressions numerically at inter- 
mediate stages and MACSYMA is able to generate fortran code. 
Collectively, these capabilities provide the means to create 
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subroutines that can be developed and verified separately. 

It is also a more efficient use of computer time to simplify 
expressions prior to numerical evaluation which MACSYMA can 
perform with variable precision on both fixed and floating 
point calculations. 

In order to illustrate an application of MACSYMA capa- 
bilities to problem solution, it is pertinent to present 
examples of manipulations related to the current study. 

The examples include the display of MACSYMA generated beam 
functions and a description of steps necessary to form the 
Rayleigh quotient for natural frequencies . The Rayleigh 
quotient is applied to a simply-supported panel utilizing 
beam sine functions as the assumed mode shapes. The example 
using sine functions was chosen for presentation because 
it is not conceptually difficult and gives an easy to follow 
MACSYMA listing. Both examples have been edited to exclude 
details of MACSYMA program statements and since MACSYMA 
exhibits no greek nomenclature, english abbreviations are 
used. The examples will also introduce the two-dimensional 
display returned by MACSYMA that is used in Appendices B, C 
and D. 

MACSYMA Ge nerate d Beam functions 

The beam functions chosen to represent the modes of 
vibration of an elastically restrained plate were given in 
equations (10) and (11). The following brief listing is a 
demonstration of computer generated expressions for these 
beam functions. The functions are inserted into equation (9) 
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and expanded in the series summation for a finite number of 
terms. 

The general expression for the assumed displacement 
functions from equation (9) is represented by MACSYMA as, 

q p 


w(x, y) 


\ 

> 

/ 


\ 

( > x(m) a (m, n)) y(n) 


n « 1 m = 1 


(27) 


The expressions for X m (x) and Y n (y) from equations (10) and 
(11) are , 


alph(m) x alph(m) x 

x (m) » b(m) sinh( ) + c(m) sin( ) 

a a 

alph(m) x alph(m) x 

+ d (m ) cosh ( — — ) + cost — ) 

a a (28) 

and , 

bet(n) y bet(n) y 

y(n) = e(n) s inh ( ) + f ( n ) sin ( ) 

b b 

bet(n) y bet(n) y 

+ g ( n ) cosh ( ) + cos ( ) 

b b (29) 

Substituting equations (28) and (29) into (27), 


\ \ 

w(x, y) * > ( > a(m, n) 

/ 7 

n * 1 m = 1 


(b(m) sinh( 


alph(m) x 


alpn(m) x 
) + c(m) sin( ) 


a 


a 
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alpn(m) x alph(m) x 

+ d(m) cosh( ) + cos ( ))) 

a a 

betiln) y bet(n) y 

(e(n) sinh( — ) + f { n ) sin ( ) 

b b 


bet(n) y bet(n) y 

+ g (n ) cosh ( ) + cost )) 

b b 


(30) 


and finally expanding to show four terms ( p=2 , q=2) we 
have, 

alpht 2 ) x 

w(x, y) - (at 2, 2) (b(2) sinh( ) 


alpht 2 ) x alpht 2 ) x 

+ c( 2 ) sint ) + d( 2 ) cosht ) 

a a 

alpht 2) x alpht 1) x 

+ cost )) + at 1 , 2) (b( 1 ) sinht ) 

a a 

alpht 1 ) x alpht 1 ) x 

+ c(l) sint ) + dtl) cosht ) 

a a 

alph(l) x bet ( 2 ) y 

+ cost-- ))) ( e ( 2 ) sinht ) 

a b 

bet (2) y bet (2) y bet (2) y 

+ ft 2) sint ) + g ( 2 ) cosht ) + cost )) 

b b b 

alpht 2 ) x alpht 2 ) x 

+ (at 2 , 1) ( b( 2 ) sinht ) + c{ 2 ) sint ) 

a a 

alpht 2 ) x alpht 2 ) x 

+ d(2) cosht ) + cost )) 

a a 
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alph(l) x alpn(l) x 

+ a ( 1 , 1) ( b( 1 ) sinh( ) + c ( 1 ) sin ( ) 

a a 

alph(l) x alph(l) x 

+ <3(1) cosh( ) + cos ( ))) 

a a 

bet(l) y bet(l) y 

( e ( 1 ) sinh( ) + f (1) sin( ) 

b b 


bet ( 1 ) y bet ( 1 ) y 

+ g(l) cosh ( ) + cost )) 

b b 


(31) 


Simply-Supported Plate Solution - 
( MACSYMA ) , Rayleigh Quotient 

The expression for the Rayleigh quotient that is formed 
from the plate strain and kinetic energies was given in 
equation (4). Another demonstration of MACSYMA capabilities 
is provided whereby the integral expressions for the energy 
terms are computer generated, the quotient is formed and 
beam functions are inserted. The partial derivatives are 
evaluated for the strain energy integrand expressions and 
integration is performed over the plate surface. 

The simply-supported plate solution can be represented 
by sine functions taken from the simply-supported beam. When 
these functions are used, the equation for natural frequen- 
cies reduces to the well known exact expression. 

Forming the quotient from equation (4), MACSYMA gives, 

b a 

/ / 2 

2 lid 2 

nfmn * a (I I (( (w(x, y ) ) ) 

radian 1 J 2 

/ / dy 

0 0 
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2 2 2 
d d d 2 

+ 2 pr ( (w(x, y))) ( (w(x, y))) + ( (w(x, y ) ) ) 


dx 


dy 


dx 


d 2 

- 2 pr ( (w(x, y ) ) ) + 2 ( 

dx dy 

b a 

/ / 

[ l 2 

/(h rho I I w (x, y) dx dy) 

J J 
/ / 

0 0 


d 2 

(w(x, y) ) ) ) dx dy) 

dx dy 


(32) 


The normal modes of vibration of a simply-supported 
plate can be represented by, 

ftpi m x %pi n y 

wmn sin ( ) s:Ln( ) 

a b 

and substituting equation (33) into (32), 


(33) 


2 

nfmn * d 

radian 

b a 

/ / 2 

[ ( d %pi m x %pi n y 2 

(I I (( (wmn sin( ) sin( ))) 

] 1 2 a b 

/ / dy 

0 0 

2 

a %pi m x %pi n y 

+ 2 pr ( (wmn sin( ) sin ( ))) 

2 a b 

dx 
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d %pi m x %pi n y 

( ( wmn sin ( ) sin( ))) 

2 a b 

dy 


d %pi m x %pi n y 2 

+ ( (wmn sin( ) sin( ))) 

2 a b 

dx 


d %pi m x %pi n y 2 

- 2 pr ( (wmn sin( ) sin( ))) 

dx dy a b 


d %pi m x %pi n y 2 

+ 2 ( (wmn sin( ) sin( ))) ) dx dy) 

dx dy a b 


a b 

/ / 

2( 2 %pi m x ( 2 %pi n y 

/(h rho wmn (I sin ( ) dx) I sin ( ) dy) 

J a ] b 

/ / 

0 0 

(34) 

The respective derivatives are taken and MACSYMA yields. 


2 

nfmn = d 

radian 

b a 422 22 %pi m x 2 %pi n y 

/ / 2 %pi m n pr wmn sin ( ) sin ( ) 

t l a b 

(I I ( 

J 1 2 2 

/ / a b 

0 0 

44 2 2 %pi m x 2 %pi n y 

%pi n wmn sin ( ) sin ( ) 

a b 


4 

b 
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44 2 2 %pi m x 2 %pi n y 

%pi m wmn sin ( ) sin ( ) 

a b 


422 22 %pi m x 2 %pi n y 

2 %pi m n pr wmn cos ( ) cos ( ) 

a b 

2 2 
a b 

422 2 2 %pi m x 2 %pi n y 

2 %pi m n wmn cos ( ) cos ( ) 

a b 

) dX fly) 

2 2 
a b 


/ 

2 ( 2 %pi m x 

/ ( h rho wmn (I sin ( 

] a 

/ 

0 

and performing the integration, 
2 


/ 

[ 2 %pi n y 

) dx) I sin ( ) dy) 

] b 

/ 

0 

(35) 


nfmn 


radian 


44 4 42322 454 

d ( %pi a bn +2 %pi a b m n + %pi b m ) 


4 5 

a b h rho 

Finally, a simplification gives. 


(36) 


nfmn 


radian 


4 2 2 2 2 2 

%pi d (a n + b m ) 

4 4 

a b h rho 


(37) 


which is the exact expression. 
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Considerations to Current Approach 

The first consideration to the Rayleigh-Ritz approach 
attempted herein is the generation of oeam functions that 
display unsymmetric boundary conditions. Although the algeb- 
raic manipulations could be done Dy hand, the intermediate 
calculations involved are long and cumbersome. Suostitution 
of the general expressions of the oeam functions into three 
of the equations of the boundary conditions involves taking 
derivatives (first, second and third order) and simultaneous 
solution for the coefficients, B, C and D. The expressions 
for these coefficients are then inserted into the fourth 
boundary equation yielding a transcendental equation for the 
argument, d • This transcendental equation is solved numer- 
ically using the Newton-Raphson method. This requires the 
derivative with respect to & of that equation. Appendix B 
displays the lengthy expressions returned by MACSYMA for 
these calculations. For comparisons of degree of difficulty, 
the calculations were performed by hand. This comparison is 
descriDed in chapter six. 

Once the beam functions are derived, it is necessary to 
obtain expressions for the integrals of the energy terms for 
natural frequencies. Again, these computations could be done 
by hand, however without symmetric simplifying assumptions, 
the manipulations are arduous. Carmichael assumed symmetric 
restraint conditions on parallel edges and was rewarded with 
concise easy to program results. By examining the necessary 
integrals with consideration paid to displacement and slope 
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on those parallel sides, the integrals were reduced to 
explicit expressions in terms of the coefficients (see refer- 
ence 13). Unfortunately, the assumption of unsymmetric 
boundary conditions in this study does not allow the reduction 
of the integrals to simple expressions. With the aid of 
MACSYMA, these integrals are directly evaluated over the 
plate surface resulting in the equations given in Appendix C. 

Carmichael demonstrated another simplification when 
considering the extraction of natural frequencies. By diagon- 
alizing the eigenvalue problem, the mn t ^ 1 mode is represented 
by the mn th term only. This was an acceptable approximation 
since the off diagonal terms are much smaller than the diagonal 
terms in the frequency determinant. As a result, the off diag- 
onal terms make minimal contribution to the calculation of 
natural frequencies. Carmichael's analysis displayed an error 
of less than 1% by using this type of approximation. Diagon- 
al ization of the frequency determinant is used in this study 
also to simplify the derivation and programming of solution. 
Accuracy of the results is shown in chapter five. 

Another consideration to the MACSYMA aided solution 
method relates to the determination of the coefficients, A mn , 
of forced response. Equation (22) shows the necessary equation 
for this calculation. It should be noted that equation (22) 
is for an undamped structure. Damping is incorporated into 
the computation in Appendix D. The damping term is inserted 
once the expression for the undamped A mn is found. The 
equation for these coefficients requires the evaluation of 
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the integral formed by the work done by the pressure on the 
plate (equation (19)). Although this integral is not as dif- 
ficult as those from tne plate strain and kinetic energies, 
MACSYMA again performs the integration conveniently. 

As a final note to the determination of the coefficients, 
a similar diagonal ization is made to the simultaneous equations 
arising from equation (22). Roussos discovered that the off 
diagonal terms are again small in comparison to the terms on 
the diagonal (see reference 5). This simplification aids in 
the programming of the solution algorithm. The accuracy resul- 
ting from this assumption is given in chapter five. 
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CHAPTER 4 

Finite Elements as a Solution Method 

* 

Introduction and Description of NASTRAN 

Since the early 1960's, the finite element method has 
had an accelerated development and has gained wide accep- 
tance to many engineering applications. Although it was 
initially designed to be used in structural analysis, it has 
been successfully employed in other engineering disciplines 
such as fluid mechanics and heat transfer. Perhaps the main 
reason the finite element method has experienced such rapid 
advancement is the computational power given by the computer. 

Essentially, the finite element method (as applied to 
structural analysis) is an energy formulation similar to the 
Rayleigh-Kitz approach. However, rather than define admissible 
functions over the entire domain, functions are chosen to 
describe the behavior over elements or subdomains. That is, 
the continuous domain is discretized into a finite number of 
elements. The functions (called element interpolation functions) 
are usually easier to define than admissible functions over 
an entire domain. In fact, low order polynomials are often 
chosen which facilitates integration of the energy terms. As 


* Acronym for NASA STRuctural AN alysis Program 


Brewer 38 


n result, the method can be applied to structures of compli- 
cated geometry. Once the equations are used to determine the 
elemental mass matrices, stiffness matrices and force vectors, 
they are assembled to describe the behavior of the entire 
structure. Compatibility of displacement is required at the 
grid points connecting the elements as the global matrices and 
vectors are assembled. The manipulations of matrix algebra 
are then used to to retrieve the solution data. 

NASTRAN is a comprehensive finite element software 
package that has been developed with general purpose objec- 
tives (17). The program can be used to examine structures 
of any size, shape and configuration. It will handle structures 
that exhibit isotropic to generally anisotropic elastic 
relations. The program is able to perform real or complex 
matrix operations and determine vibration frequencies and 
modes. Various loadings may be applied to the structural grid 
points including concentrated loads and distributed loads. 

The loads can be transient, sinusoidal steady state and 
random. Flexibility nas oeen incorporated into the NASTRAN 
program package in an effort to anticipate changing needs 
and applications. As described, NASTRAN is quite capable to 
analyze uniform rectangular plates with elastic restraining 
springs on the edges. Rectangular plates are conveniently 
modeled with a cartesian coordinate grid configuration. 

NASTRAN uses the displacement approach with reference to the 
grid points to analyze these structures. NASTRAN uses a 
polynomial representation of the displacements. 
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A typical NASTRAN analysis involves an input file that 
is composed of three parts (18). These are the executive cont- 
rol, the case control and the bulk data. The executive control 
portion defines the type of solution (i.e. static, dynamic, 
heat transfer, etc.) to be used and estimates the time of 
computation. Several options for special features such as 
diagnostic requests are also designated in the executive 
control . 

The case control provides more detail of the input to 
the specific problem under consideration. The selection of 
eigenvalue extraction method, the specification of loading 
cases to be examined and the output frequency selections are 
made in the case control. Model and material properties such 
as damping of the structure are referenced and grid point cons- 
traints are selected from the bulk data. Requests of type, 
location and sorting of the output are issued in this section. 

The third and most detailed part of the NASTRAN input is 
the bulk data. This section is a formatted data list that des- 
cribes the model and material properties. The grid configuration 
is specified by grid point location and element connectivities. 
Loading is indicated with amplitude and time or frequency 
variation as well as points of load application. Parameters 
that control accuracy of solution are listed. The bulk data 
is the "meat" of the NASTRAN input and can be a long listing. 

As a result, it is the bulk data that presents the best chance 


for format error. 
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Figure (2). NASTRAN Finite Element Model 
Description of Elements and Models 

The models (figure (2)) chosen for analysis in this study 
were composed of isotropic quadrilateral (CQUAD1) elements of 
uniform thickness. The CQUAD1 element possesses both in-plane 
and bending stiffnesses although the in-plane motions were 
neglected since the two motions are uncoupled in small deflec- 
tion linear transverse vibrations. This element uses two sets 
of coplanar overlapping bending triangles (figure (3)) that 




Figure (3). CQUAD1 Element Consisting of Overlapping 
Pairs of Triangular Plate Elements 





Brewer 41 


each have one half the quadrilateral thickness. Line segments 
joining points are assumed rigid. The stresses are computed 
for each triangle at the intersection of the diagonals and 
averaged (19). For rectangul-ar elements, this calculation is 
at the center of the element. Strains are not calculated in 
NASTRAN dynamic analysis. The accuracy returned by the CQUAD1 
element is described in reference (19). 

The other element necessary is the scalar spring element 
to model elastic boundary restraints. Scalar spring elements 
( CELASl ) were attached to each grid point on the four edges 
of the plates. Two spring elements were fixed to each grid 
point (one for translation and one rotation). At the corners, 
however, it was necessary to attach an additional rotational 
spring in order to model the joining of two respective edges. 
The connections were completed by fixing each spring from 
grid point to grounded scalar points. The bulk data was 
arranged so as to allow the assignment of the eight distinct 
spring constants (four edges with two types of spring on each 
edge ) . 

Six models were employed to conduct the study. They 
included 15"x 3", 15"x 6", 15"x 9", 15"x 12" and two square 
plates at 12"x 12" and 15"xl5". These models gave aspect 
ratios (=b/a) of 0.2, 0.4, 0.6, 0.8 and 1.0, respectively. 

The input data given to NASTRAN for each model was represen- 
tative of the material properties of a typical aluminum 
panel. Upon retrieval of the output data, the frequencies 
and responses were nondimens ionali zed as described in the 
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results. The analysis was conducted using NASTRAN normal 
modes frequency response. 

Eiqenva 1 ue Extraction 

Once the global mass and stiffness matrices are assem- 
bled, the eigenvalues and eigenvectors are extracted accor- 
ding to the formula, 

( [K] - X CM] ) { <£} =0 (38) 

where, [M] and [K] are the global mass and stiffness matrices, 
X and I (p } are the desired eigenvalues and eigenvectors, 
respectively. The method selected for this extraction was 
the Inverse Method with Shifts (19). The NASTRAN software 
package contains the algorithm for this procedure. 

The inverse method is a vector iteration procedure. 
Writing equation (38) as, 

W} n+1 = [D] { u ) n (39) 

where , 

Id] = [ K ] ~ 1 [M] (40) 

and lul n is the trail input vector (arbitrary initial trail 
vector, |u}|) and { v ] n+ ^ is the vector obtained from the 
n t * 1 iteration. 

Then the n+] st trail vector input to tne iteration is, 

< u >n+l = X(n+1) < v >n+l < 41) 

where X(n+1) an a PP r °P r i ate scaling factor. NASTRAN 

selects X(n+1) as t ^ e i- nverse °f the element of largest 
absolute value in (v} n+1 so as to normalize {u} n+ i with +1 
as the largest element. It has been proven (19) that (v} n+ j 
converges to the eigenvector, { <p > and X( n+ u converges on 
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the corresponding c*i genva 1 ue , \ . In most cases, the procedure 
results in the lowest eigenvalue. 

In order to obtain eigenvalues and eigenvectors of higher 
modes, the algorithm is employed with spectrum shift. That 
is, a shift point, X C) is selected that is close to a higher 
eigenvalue and a shifted eigenvalue is defined as, 

A - X- X 0 l42 > 

Substituting into equation (38), 

< [ K* ] - A f M] ) { } =0 (43) 

where, 

t K* ] = [K] - X Q [M] (44) 

The inverse method is then applied to this shifted 
eigenvalue problem. A is solved for and the required eigen- 
value is, 

X - Xo + A (45) 

The convergence criteria used by NASTRAN is based on ortho- 
gonality. A check of retrieved eigenvectors is made such that 
the modal mass matrix, 

{ Cpi 1 T [M] { <£j } <7 i = j (46) 

where Y is a user supplied parameter and in the ideal case, 

r = o. 

Finally, mention should be made about the use of incon- 
sistent (lumped) and consistent mass formulations. The lumped 
mass formulation assembles mass matrices by lumping equal 
amounts of mass at the grid points. As a result, the mass 
matrix is diagonal. This type of simplification tends to 
lower the extracted eigenvalues similar to the effect of a 
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larger denominator in the Rayleigh quotient (equation (4)). 
Although this may seem favorable- by the argument of upper 
bounded eigenvalues, computational accuracy may suffer. The 
natural frequencies obtained in this formulation converge on 
the true values from below as the number of elements increases. 

In the consistent mass formulation, the off diagonal 
terms are nonzero. These terms couple the adjacent grid 
points. The eigenvalues retrieved in this case converge more 
rapidly with a smaller number of elements. 

Inadvertently, initial program runs using NASTRAN were 
done under the lumped mass conditions. After subsequent 
checks of the results from NASTRAN employing the consistent 
mass formulation, it was found that the finite element models 
had been constructed with a sufficient number of elements 
for good convergence of natural frequencies of lower modes. 

Al 1 NASTRAN runs were done using consistent mass thereafter. 
Load i n g and G eometri c Considerations 

To properly distribute the continuous uniform elastic 
boundary restraints along each edge on the finite element 
mode 1 s , equivalent discrete values must be determined. A 
scalar spring was attached to each edge grid point and values 
were assigned according to the formula, 

R, or K i mastran = ( R i or K i / unit length) x side length 

I elements along side 

(47) 

figure (4) shows how the values were concentrated at a typical 
edge grid point. Special consideration was given to the 
corner points for each respective side which necessitated 
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the assignment of values equal to one half those from 


equation ( 47 ) . 

^ K a or R a ( 


Continuous 


'//////////////A 


or R^ Discretized 




Length of Continuous Spring Length Covered by Discrete 

„ / „ % Spring Value 

Figure (4). * 

Similar to the discretized spring values, the uniform 
pressure load must be divided and concentrated for normal 
incidence application at every grid point. The equivalent 
grid point load amplitude values were determined as, 

P NASTRAN * i-S / 7 t ° t ; 1 area 

total # of elements 


Figure (5) displays the equivalent geometric application of 


Pressure, P Continuous 
(Uniform) 


77 // 77 / 
////// / 
/////// 
/////// 

Y//AYA 


Area of Uniform Loading 


Pressure, P Discretized 



ig Area Covered by Discrete 

Loading Values 

Figure ( 5 ) . 
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load at typical grid points. Again, special consideration is 
required. The edge grid points carried a load amplitude of 
one half the values from equation (48) and the corner points 
were assigned 1/4 values. 

The accuracy of results returned by this finite element 
model was previously demonstrated (5). 
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CHAPTER 5 

Results and Comparisons 

Frequency and Restraint Parameters 

In order to assess the validity of tne single term Ritz 
solution derived using elastically supported beam functions 
and the results given by NASTRAN , comparisons are now made to 
frequency parameter studies presented in the literature under 
various selected edge conditions and aspect ratios. There is 
a vast amount of literature available on the determination of 
frequency for the classical simply-supported, clamped and free 
boundary conditions and the 21 possible combinations from the 
mixing of these three among the four edges. In general, those 
cases are considered special limiting cases for the current 
study. Mainly, comparisons are made to the published results 
that address the problem of elastically restrained plates as 
most have been shown to converge on accepted values of fre- 
quency for the classical cases. Therefore, the classical edge 

★ 

conditions are noted only as they occur in the limiting case. 
The available literature on vibration of plates with elastic 
restraints is sparse, however, especially for plates possessing 
both translational and rotational springs at the boundaries. 

* Numerical values for zero and infinite spring stiffness 
used in this study were 10~ 4 and 10 + , respectively. 


Brewer 


48 


So that the results may be generalized and considered 
for all uniform isotropic panels, the frequency and restraint 
parameters are non-dimensional ized according to the formulas 
given by Warburton and Edney (12). The dimensionless frequency 
parameter is defined as, 

= U) b 2 Jpn7o (49) 

and the rotational and translational restraint parameters 
are given as. 


* 

R = Rb/D x=o , a 

= Ra/D y=o , b 

* 7 

K = Kb^/D x=o , a 

= Ka 3 /D y=o , b 


(50) 

(51) 


The method given by Carmichael (13) is generally regarded 
as an acceptable Ritz approximation to determine natural freq- 
uencies for the cases when edges are between simply-supported 
and clamped (ie. K=inf as R varies). Table (1) shows the com- 
parison of frequency for selected values of restraint. The 
first three modes at three values of restraint parameter for 
five aspect ratios are displayed. It is seen that the new 
elastic Ritz approximation tends to default to Carmichael's 
results to within 1%. The discrepancy between the finite 
element solution and the results of both Ritz solutions is 
slightly greater and tends to increase as the aspect ratio 
decreases. The maximum percent difference is 7% occurring at 
an aspect ratio of 0.2. 

Figure (6) shows the variation of fundamental frequency 
parameter for a square plate as both types of restraint vary 
from values of zero and infinity as given by Warburton and 
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Comparison of Characteristic Frequency Mran**t«T 

inf , Rj *l<2 * Kj*U4 *R 
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inf 
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19.39 

f 20.17 

21.54 
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20 

19.39 
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19.71 
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22.66 

23.49 
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inf 

22.59 

23.42 

24.85 

* * 



23.14 

24.86 

25.62 

• ** 


• Carmichael (Kef. 13) 
** Hits Solution herein 
*** Nastran F.E. models 


Table (1) 

Edney. All boundaries are subject to equal restraint. The 
curve labeled R* =S and K* =inf shows variation from simply- 
supported to clamped edges and the lowest curve marked R* =0 
and K =S displays the changes between the limiting free and 


simply-supported boundary conditions. Intermediate to these 


two curves (marked R* ~S and K* =S) is the change in frequency 
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as the edges of the plate vary from free to clamped conditions 
as both rotational and transverse boundary stiffnesses are 
increased simultaneously. As stated previously, the assumed 
displacement functions used by Warburton and Edney consisted 
of a weighted summation of beam functions displaying the 
simply-supported, clamped, free and sliding end conditions. 



Figure It). Variation of fundamental frequency parameter 
JXjj of a square plate with translational and rotational 
restraint parameters. Identical boundary conditions on 
all edges. (Taken from reference (121). 

In comparison, figure (7) displays the same study of 
variation of fundamental frequency as returned by the Ritz 
solution assuming beam displacement functions that satisfy 
the general elastic edge conditions and results given by the 
NASTRAN finite element solution. The top and bottom curves, 
corresponding to changes in one type of boundary spring while 
the other is held constant, show good agreement to the results 
of figure (6). Indeed, the results of all upper and lower 
curves agree to within 1% for both figures. The intermediate 
curve of figure (7), however, displays a notable discrepancy 
when compared to figure (6). The elastic Ritz and finite 
element solutions agree to within 5% but both return decidedly 
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lower values of frequency for all values of S on the middle 
curve except for the limiting cases of the free and clamped 
boundary conditions. 



Piaur* (7). Variation of fundamental frequency parameter 
of a square plate with translational and rotational 
restraint parameters . Identical boundary conditions on 
all edges. 

Returning to discrete numerical comparisons, Filipich, 
Reyes and Rossi (6) tabulated various fundamental frequency 
results for selected values of rotational and translational 


restraint. They assumed functions given by, 


W(x,y)= 



(52) 


where, 

x i = Cli x 4 /a 4 + x 2 /a 2 + 1 (53) 

and, ' 

y i = costn^ y/b) (54) 

and the computations were executed using the Galerkin method. 

Table (2) shows these results for the case of equal rotational 

restraint on all edges for a square plate while the transverse 

edge restraint is held to infinity. This case corresponds to 

the top curves of figures (6) and (7). For further comparison. 
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the results given by Mukhopadhyay (7) who used a semi-analytic 
solution method by applying rotational ly restrained beam 
functions in one coordinate direction and a finite difference 
approximation. in the other direction. The limiting simply- 
supported and clamped frequency parameter values from Leissa 
(11) are also included. Table (2) shows general agreement 
between all solutions, however a maximum percent difference 
of 10.5% exists between the results of references (6) and (7). 
The elastic Ritz solution and the finite element results 
differed from all by no more than 5%. 


Comparison of fundamental frequency parameter 
Square plate 

*l“ K 2" K 3" K 4“ in * Rj"R2*H 3 »H 4 »R 


Rb/D 

inf . 

100 

10 

1 

0.1 
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_j 
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19.74 
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29.55 

21.51 

19.94 

19.74 

1 

35.76 

34.46 

29.32 

21.32 

19.75 

19.56 

ED 


* Filipich, Reyes a Rossi (Ref. (6)) 

*• Mukhopadhyay (Ref. (7)) 

Ritz Solution herein 
••** Nastran F.E. Solution 
Leissa (Ref. (ID) 

Table (2) 


Table (3) gives results that correspond to the curve 
marked R* =0 and K* =S on figures (6) and (7). Although there 
was good agreement demonstrated on the lower curves of those 
graphs, the numerical data of table (3) yields a discrepancy 
of as much as 50% between the Galerkin method of reference 
(6) and the two solutions used in this study. Oddly, Warburton 
and Edney only mentioned the comparative results stating 
qualitatively that there was good agreement between their 
results and those of Filipich, Reyes and Rossi. No numerical 
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data was given to substantiate the statement. Again, the 
agreement between the elastic Ritz and finite element solutions 
for this case is to within 1% and the closeness of these 
results to those of Warburton and Edney for R* =0 and K* =S 
is seen on the bottom curve of figures (6) and (7). 


Comparison of fundamental frequency parameter 
Square plate 

Rj*R 2 *R3*R4»0 K l* K 2* K 3 mK 4* K 


Kb 3 /D 

inf 

11000 
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10 

1 

fl 

HI 


19.74 

IQS 
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ii. je 

5.92 

n 

HI 

XL 1 1 
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Rlu 
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5.97 
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IS 


19.56 
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14.44 

5.99 

1.99 

■fl 

mm 


• Filipich, Reyes t Rossi (Ref. (6)) 
** Ritz Solution herein 
**• Nastran F.E. Solution 


Table (3) 


As a final example of equal restraint on all of the four 
edges, the rotational restraining springs are held at infinite 
stiffness rather than zero shown in table (4). The transverse 
edge restraint is allowed to vary and the fundamental freq- 
uency parameter is tabulated for selected values. This case 
may be considered as a sliding boundary on foundation springs. 
The results display a similar disagreement of frequency between 
the Galerkin solution and the two solutions of this treatise. 
Filipich, Reyes and Rossi achieved frequency calculations 

that were higher than the elastic Ritz and finite element 

* 

results with a 60% discrepancy at K =10. As before, the elastic 
Ritz and NASTRAN results were in close agreement with no more 
than 5% difference noted. 

Tables (5), (6) and (7) address the cases where parallel 
sides of a square plate are restrained to the same degree. The 
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Comparison of fundamental frequency parameter 
Square plate 

B l“ B 2* B 3“ B 4* : *’ nf K l -K 2 «K3»K 4 »K 
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• Filipich, Reyes t Rossi (Ref. (6)) 

** Ritz Solution herein 
*** Nastran F.E. Solution 

Table (4) 

solutions used for comparison and the variations of boundary 
spring parameter are similar as for tables (2), (3) and (4), 
respectively. First, the transverse edge restraint is held to 
infinite stiffness as changes are made to the rotational 
restraint. The variations of fundamental frequency are disp- 
layed in table (5). The results of the top row correspond to 
the case where one pair of opposite edges are clamped while 
the rotational restraint on the remaining parallel sides is 
relaxed to converge on the case of opposite sides clamped and 
the remaining boundaries are simply-supported. An exact value 
results when parallel edges are simply-supported and is taken 
from reference (11) for comparison. The bottom row of table 
(5) shows the variation of frequency as two opposite boundaries 
are held at the R=0 (simply-supported) value and the remaining 
two sides change from the clamped (R=inf) to. the simply- 
supported condition. All solutions show good agreement and 
display at most 3% difference. 

Presented in table (6) is the variation of frequency as 
parallel edges are changed from simply-supported to free end 
conditions. The rotational restraint is set to zero and the 
transverse boundary spring is allowed to vary numerically 
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Comparison of fundamental frequency parameter 
Square plate 
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• Filipich, Reyes i Rossi (Ref. ( 6 )) 
Mukhopadnyay (Ref. (7)) 

*** Rits Solution herein 
Nastran F.E. Solution 
Leissa (Ref. (Ill) 

* Table (5) 


from infinity to zero. The elastic Ritz solution shows a more 
rapid descent than either the Galerkin solution of reference 
(6) or the NASTRAN frequency predictions. Intermediate to 
the limiting cases of all edges simply-supported or free is 
the case of opposite sides free while the other parallel 


edges are simply-supported. This case has an exact value and 


is given by Leissa (11) as 9.6314. Although the values of the 


elastic Ritz and Galerkin solutions tend to converge on this, 
the finite element solution gives an inflated value with a 
percent difference of 76%. Convergence of the NASTRAN results 
is exhibited for only cases of equal restraint on all edges 
in table ( 6 ) . 

For the study shown in table (7), the rotation on the 
sides is completely constrained (ie. R=inf) while the trans- 
lational boundary restraint is changed for equal values on 
parallel edges between infinity and zero. The results show 
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good agreement between the elastic Ritz and finite element 
solutions, however both are much lower than the results of 
reference { 6 ) . 


Comparison of fundamental frequency parameter 
a Square plate 
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Comparison of fundamental frequency parameter 
Square plate 
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Table (7) 

Tables (8), (9) and (10) display comparisons of charac- 
teristic frequency parameter for selected cases given by 
Mukhopadhyay when edges are restrained to different degrees. 
The boundaries of the plate are completely constrained in the 
transverse direction ( ie . K=inf). As such, edge conditions 
lie between the simply-supported and clamped states. Each of 
the tables presents results for plates with boundary springs 
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on only one side at x=a (ie. side marked 2 on figure (1)). 

The remaining three sides are held at either the clamped or 
simply-supported condition. Frequency results are given for 
four modes at three values of restraint parameter. The aspect 
ratios selected are for b/a=1.0, 2.5 and 0.4. The ratio of 
2.5 was chosen since, for the cases of unsymmetric boundary 
restraints, natural frequency is dependent on the magnitude 
of restraint and the side length upon which it is applied. 
This ratio was achieved by interchanging the x and y coor- 
dinate axes on the plate with b/a=0.4. 

In table (8), three sides of the panel are held at the 
simply-supported state while the fourth edge is subject to 
elastic rotational support. Three select values of restraint 
parameter are displayed. The agreement of comparison of 
results of reference (7) are good, in general. The NASTRAN 
finite element solution gives slightly lower values for all 
cases other than the first mode at b/a=0.4. The elastic Ritz 
solution agrees more closely with the results of Mukhopadhyay . 
The maximum percent difference noted is 3%. 

For the cases presented in table (9), three of the sides 
are clamped as the fourth edge is restrained elastically in 
rotation. Again, the results of all solutions show general 
agreement, however for higher modes at b/a=2.5, the elastic 
Ritz and finite element frequencies are higher than those of 
reference (7). The maximum percent difference is found there 
and is approximately 4%. 
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Comparison of Characteristic Frequency Parameter 
K i *K 2 *K3*K 4 *mf R 1 3**4*° k2=R 


b/a 

Rb/D 

mode 1 

mode 2 

mode 3 

mode 4 


1.0 

0.001 

19.46 

19.74 

19.56 

49.21 

49.35 

48.90 

49.35 
49. 36 
48.91 

79.88 

78.95 

77.65 

• 

* # 
** » 

1.0 

20.19 

20.19 

20.00 

49.32 

49.54 

49.08 

50. 15 
50.09 
49.67 

~ 19773 ” 
79.42 
78.12 

”1 

• * 

* * * 

lOOoO 

23.60 

23.44 

23.17 

5T73Q 

51.77 

50.85 

5TTF“ 

57.75 

57.57 

" 4578 9 

85.39 

84.91 

• 

* • 
**• 

2.5« 

0.001 

71.66 

71.55 

71.51 

101.33 

101.16 

100.02 

150.04 

150.49 

147.12 

219.57 

219.57 

214.18 

• 

* • 

* * * 

1.0 

76.30 

73.57 

73.55 

104.73 

102.60 

101.49 

152.41 

151.48 

148.14 

218.36 

220.25 

214.90 

* 

* * 

• • * 

100.0 

"1 02.31 
99.73 
100.15 

127.08 

124.48 

124.01 

1 70. 99 : 

169.24 

166.32 

233.99 

235.01 

229.42 

* 

** 

0.4 

0.001 

11.39 

11.45 

11.44 

16.16 

16.19 

16.01 

24.11' 

24.08 

23.54 

35.24 

35.13 

34.27 

• 

• * 

• * * 

1.0 

11.42 

11.51 

11.49 

16.22 

16.32 

16.13 

24.21 

24.29 

23.75 

35.35 

35.39 

34.56 

• 

*« 

»** 

100.0 

11.67 

11.81 

11.71 

17.06 

17.17 

16.91 

25.70 

25.81 

25.25 

37.51 

37.62 

36.95 

• 

** 
*• * 


• Mukhopadhyay (Ref. (7)) 

Rit* Solution herein 
**• Nastran F.E. Solution 

Table (8) 

Table (10) gives results for the case of a plate pos- 
sessing clamped boundary conditions on two adjacent edges. 


one of the remaining edges is simply-supported and the last 


is elastically restrained in rotation. The results show 
similar agreement as tables (8) and (9). It is noted that 
the elastic Ritz and NASTRAN solutions give frequencies that 
are not consistently higher or lower but may be either. The 
maximum percent difference is approximately 6%. 

The final comparison of frequency is given in table (11) 


Since no data is available in literature for unsymmetric 
elastic boundary conditions with both translational and 
rotational springs, tnese results compare only the elastic 
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Ritz and NASTRAN finite element solutions. The fundamental 
frequency parameter is tabulated for selected values of 
restraint for a square plate. The rotational and translationa 
dimensionless restraint parameters are held equal in mag- 
nitude at 200. Edges are clamped unless specified with 
restraining values. Except for the case of all edges clamped, 
the finite element frequency prediction is higher than the 
elastic Ritz results. The agreement' of the two is generally 
good, however an 8.4% difference is seen for the case of 
parallel edges elastically restrained with the remaining two 


boundaries clamped. 


Coaparison of Characteristic Frequency Parameter 
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Comparison of Characteristic Frequency Parameter 
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Displacement , Acceleration and Strain 

The response calculations necessary for design and control 
of flexible structures is not limited to frequency prediction. 
Although discrepancies exist between all methods used to 
predict frequency and academic agreement on an accepted num- 
erical or analytical approach and results can take years, the 
practicing engineer needs prediction methods for measurable 
physical response such as displacement, acceleration and 
strain in order to compare with testing and real applications. 
Of course it is unreasonable to compare physical response 
given by a method under any dynamic load if the frequency 
predictions are unrealistic or disagree. To this end, sample 
comparisons of the physical responses are now given. 

It was stated earlier tnat the loading is assumed to be 
a series of harmonic inputs to the mechanical system. The 
additional assumption of linear small deflection plate theory 
allows the use of superposition. That is, if the loading can 
be represented by a harmonic series, the response from each 
term in the series may be calculated individually and summed 
over the series for total response. This summation can also 
be carried out separately for contributions to total response 
from otner modes of vibration that may be excited. However, 
the fundamental mode is the easiest to excite and fundamental 
response makes the major contribution to overall response. 

This is not to say that the contribution from higher modes 
should not be considered, but in many cases, it can safely 
be neglected. For the results presented, loading is taken as 
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a spatially uniform discrete sinusoidal pressure with a 
frequency equal to the fundamental frequency of the plate. 

The single mode harmonic steady-state responses for disp- 
lacement, acceleration and strain are given. Each of these 
responses is inversely proportional to structural damping. 
That is, if damping is doubled, the response is cut in half. 
The damping used in this study is 5% critical damping. 

The solution of steady state harmonic response of a 
rectangular plate subject to a uniform pressure loading given 
by the MACSYMA aided elastic Ritz formulation is continuous. 
The root mean square displacement, acceleration and strain 
expressions represent a continuous distribution over the 
plate domain. By specifying coordinates of interest, these 
response quantities can be examined anywhere on the plate. 

As well, the strain may be calculated' through the panel 
thickness with the specification of the z coordinate. The 
assumption of linear plate theory and the absence of shear 
deformation give a linear distribution of strain through the 
thickness. Maximum values of strain are found at the top and 
bottom surfaces of the plate since linear theory assumes that 
the midplane of the panel is the neutral surface of zero 
stress/strain. The top and bottom surface strain values are 
of equal magnitude. 

The finite element formulation of NASTRAN does not give 
a continuous distribution of response quantities. Since the 
domain is divided into discrete subdomains, the distribution 
of these numerical quantities is discontinuous. The 
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displacement and acceleration responses are calculated at 
nodal points at the corners joining elements while the stresses 
are obtained from the center of the rectangular quadrilateral 
elements. These stresses are assumed constant over the element. 
The linear distribution of stress through the thickness of 
the plate is retained by NASTRAN, however. 

In order to make fair comparisons of response between 
the elastic Ritz and NASTRAN results, the plate domain was 
discretized on the finite element models using a cartesian 
mesh configuration. This configuration places a node point 
exactly at the center of the plate enabling direct comparison 
of center displacement and acceleration. Without mesh refine- 
ment, center stress is not directly calculated by NASTRAN. 
Instead, comparisons are made at a point in the center of an 
element nearest the plate center. The strain calculations 
done using the elastic Ritz formulation are adjusted to 
coincide with the coordinates of this near center point on 
the models. Reference to figure (2) on page 39 clarifies the 
location of these points. 

It was stated earlier that NASTRAN dynamic analysis 
does not calculate strain. Stresses are returned from cal- 
culations requested on a given element. For comparison, these 
stresses are converted to strains assuming plane stress by 
the formulas, 

€ X = ( 0*x- fLCTy)/E (55) 

€y = ( (T y- /X(T X )/E 


( 56 ) 
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Upon retrieval of the responses from both methods of 
solution, the displacement, acceleration and strain values 
are non-dimens ionali zed. This is done for general comparison 
to all isotropic uniform rectangular plates subject to a 
uniform sinusoidal pressure load at the fundamental frequency. 
The dimensionless displacement is given by, 

w* = (D w c )/(P b 4 ) (57) 

and the dimensionless acceleration is, 

a* = (D a c h)/(P b 4 g) (58) 

and finally the null -dimensional strain can be given by, 

€* = (2 D (l+^l) 6 c )/(h P b 2 ) (59) 

Figures (8) through (13) may best be interpreted by 
referring to figure (7). The comparison of central pnysical 
responses is given in the following figures for fundamental 
resonance excitation of a square plate with identical boun- 
dary conditions on all edges. These graphs may be considered 
crossplots of the displacement, acceleration and strain as 
they relate first to the boundary restraint values and second 
to the corresponding changes in fundamental frequency. The 
three zones as given in figure (7) are also marked on figures 

(8) through (13). In order to distinguish between the zones 

* * * * 

marked R =S K =inf and R =0 K =S, the transition point is 
noted as the frequency of the simply-supported frequency, Xl ss * 
Figures (9) and (12) were separated into individual graphs 
according to the zones for clarity. 

The results displayed in figures (8) and (9) show com- 
parison of center displacement between tne single mode elastic 
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Ritz and the NASTRAN finite element solutions. First, figure 
(8) gives the variation of displacement with respect to res- 
traint parameter. As expected, the center displacement inc- 
reases as the elastic restraints are relaxed. In figure (9), 
the displacement is crossplotted against the range of funda- 
mental excitation frequencies encountered from the restraint 
relaxation. These two figures taken together show how center 
displacement varies over the entire domain given in figure 
(7). There is very good agreement between the elastic Ritz 

and NASTRAN finite element solutions as seen on the graphs. 

* l 

The results have been truncated at w c =10 A since these ap- 
proach displacement values outside the region that is nor- 
mally considered as small deflection. 



Figure (I). Variation of cantor displacement with rospact to 
restraint parameter, S, of a square plate with translational 
and rotational restraint parameters. Identical boundary 
conditiona on all edges. Sinusoidal load at fundamental 
frequency, lljj. 


Figures (10) and (11) show the same type of crossplots 
for center acceleration. Good agreement between the NASTRAN 
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rigura (9). Variation of contar displacement with respect 
to fundamental frequency parameter, jtn* ° * ■ square plate 
with translational and rotational restraint parameters. 
Identical boundary conditions on all edges. Sinusoidal load 
at fundamental frequency, f^j. 


and elastic Ritz results is again seen and may be expected 
since the frequency and displacement predictions were in 
close agreement. The magnitude of acceleration can be cal- 
culated by multiplying radian frequency squared times the 
displacement magnitude. Therefore, acceleration calculations 
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reflect a type of combined error for frequency and displace- 
ment computations. 



Figure (10). Variation of canter acceleration with respect 
to restraint parameter, S, of a square plate with trans~ 
lational and rotational restraint parameters. Identical 
boundary conditions on all edges. Sinusoidal load at 
fundamental frequency, fjj. 



Figure (11). Variation of center acceleration with respect 
to fundamental frequency parameter. XI 11* °* * square plate 
with translational and rotational restraint pararet«rs. 
Identical boundary conditions on all edges. Sinusoidal load 
at fundamental frequency, fjj. 


Figures (12) and (13) show the variation of strain at 
a point near the center of the plate. Figure (12A) displays 
the variation as the plate boundaries move from the clamped 
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(12A) 




figure (12). Variation of eenter strain with respect to 
restraint parameter, S, of a square plate with trar.slationsl 
and rotational restraint parameters. Identical boundary 
conditions on all edqes. Sinusoidal load at fundamental 
frequency. f|j. 

to the simply-supported states. The agreement in that region 
is good. In figures (12B) and (120 it is seen that as the 
plate edges are relaxed toward the free state, the agreement 
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between the two solutions deteriorates. As the elastic res- 
traints are relaxed to values below 10^, the NASTRAN solution 
gives decidedly lower results. There is , however, agreement 
in trend. That is, the magnitude of strain tends to rise and 
fall in similar fashion for both solutions. This agreement 
is more apparent in figure (13). The shapes of the closed 
curves are very similar in the prediction of high and low 
values of strain. 



Fiqur* (13). Variation of canter strain with raapact to 
fundamental frequency parameter, Hu* °* • •*««« pl«te 
with translational and rotational restraint parameters. 
Identical boundary conditions on all edqes. Sinusoidal load 
at fundamental frequency, fjj. 


The final comparison of response is taken from the res- 
ults given in reference (5). In that work, the search for 
causes of a large bias error between experimentally measured 
and analytically predicted strain response led to the exami- 
nation of the effect of both types of boundary springs on the 
response. For that study, only the finite element method was 
used to model a plate whose edges were elastically restrained 
in rotation aiid transverse translation since no other method 
was available at the time. The response was calculated for 
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various combinations of selected values of restraint while 
maintaining a constant fundamental frequency. At one end of 
the scale, the plate possesses infinite rotational boundary 
spring stiffness with a finite value of transverse boundary 
spring stiffness. At the other end of the scale, the reverse 
is true. Combinations of finite values for both types of 
boundary spring fill in table (12). 

The NASTRAN finite element results as given in reference 
(5) have been non-dimensionalized and the elastic Ritz sol- 
ution has been used for comparison. The fundamental frequency 
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prediction for given values of restraint varied by no more 
than 5%. The results of table (12) show good comparison for 
center displacement which seems to remain relatively constant. 
The near center strain given by NASTRAN tends to drop slightly 
while the elastic Ritz solution gives relatively constant 
values. The near edge strain from a point midway along one 
of the edges increases as the translational restraint relaxes 
and the rotational restraint is tightened. This comparison 
validates the conclusions of reference (5), since the intro- 
duction of transverse motion at the boundaries could not 
account for the bias error. 
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CHAPTER 6 
Conclusions 

Conclusions from Comparisons 

Many authors have addressed the problem of free vibration 
of uniform rectangular plates. An informative survey of recent 
studies was given by Leissa (20). Most of the studies were 
restricted to plates possessing the classical free, clamped 
and simply-supported edge conditions and in some cases, the 
academic agreement on accepted theoretical frequencies is 
still open to debate. 

Many methods have been derived and used in an effort to 
calculate more accurate frequencies with ease of computation. 
With the ability of the computer to perform complicated cal- 
culations at high speed, the ease of computation is not as 
important as it once was. In fact, one of the reasons early 
derivations of solution were simplified was to yield expres- 
sions that were easier to use in the absence of efficient 
computer methods. The worx given by Carmichael (13) in 1959 
is a case in point. The assumption of symmetric boundary 
conditions with equal rotational restraint on parallel edges 
gave a simple explicit expressions for difficult integrals 
an tne application of the Rayleigh-Ritz method. 
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The Ritz method is one of the more popular approximations 
used in plate vibrations since it has proven to be a straight 
forward valid approach. The accuracy returned by this method 
depends in large part on the functions chosen to represent 
displacement. In general, the accuracy of results can not be 
estimated with certainty. The frequencies can only be said 
to be higher than the true values. By satisfying the shear 
and moment boundary conditions, when applicable, the rate of 
convergence is aided. For plates without exact solutions (ie. 
no parallel simply-supported edges), the convergence rate 
is further enhanced by the use of beam functions. The indi- 
cations are that the "best" estimates of frequency result 
when beam functions are selected so as to satisfy the edge 
conditions exactly. This choice of functions renders good 
approximations of frequency when only the single term (dia- 
gonal terms) expressions are used. Even frequencies of 
higher modes can be estimated with suitable accuracy using 
the single term solution. 

For example. Young (21) discussed the use of the Ritz 
method for plate vibrations and assumed beam functions to 
represent displacements. For the square cantilever plate which 
is a case of unsymmetric boundary conditions, he chose the 
free-free and free-clamped beam functions for the two coor- 
dinate directions, respectively. Frequency calculations were 
Dased on an 18 term series. For the first mode. Young gave 
the value of the dimensionless frequency parameter as 3.494. 

A single term approximation based on the same beam functions 
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(jives the fundamental parameter as 3.515 for a percent dif- 
ference of 0.6%. On the fourth mode, the 18 term solution 
yields a value of 27.46 and the single term approximation 
gives 28.78 for a 4.8% difference. Nassar (22) extended the 
single term approach to include plates with rotationally 
restrained edges. The boundaries were assumed to be unsym- 
metric in restraining value. Functions were chosen as a 
weighted sum of clamped and simply-supported beam functions. 

A single term solution was also generated assuming beam 
functions of an unsymmetric rotationally restrained beam. 

Both solutions gave fundamental frequency predictions that 
varied from published results by no more than 0.5%. 

In this study, a single term solution is again applied. 
The plate was assumed to be elastically restrained in both 
rotation and transverse translation. The edge spring values 
were assumed unsymmetric. For the comparisons of character- 
istic frequency parameter, it was seen that good results 
were achieved for all cases when the edge’ translation was 
completely constrained (ie. K=inf ) . In these cases, the edge 
conditions are between the simply-supported and clamped 
states, inclusively. When compared to the solution given by 
Carmichael in table (1), both the new elastic Ritz and NASTRAN 
finite element solutions gave good results. If an exact value 
of infinity for translational edge stiffness were applied to 
the new elastic Ritz solution, it would become equivalent to 
Carmichael's approximation. 
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The further comparisons of frequency parameter for plates 
rigid in edge translation showed acceptable agreement to the 
published results. The elastic Ritz and finite element results 
were equally valid for cases of unequal rotational restraint. 
For these cases, presented in tables (2), (5), (8), (9) and 

(10), the predictions of this study never differed by more 
than 1%. Tne frequency parameter comparisons made to available 
data for higher modes were within 5%. This indicates that the 
new results were soundly derived from valid approaches. 

When both types of spring were activated simultaneously 
at the plate boundaries so that transverse edge motion was 
introduced, discrepancies in frequency prediction were encoun- 
tered. Referring to figures (6) and (7), it was shown that 
the elastic Ritz and NASTRAN solutions agreed with the results 
of Warburton and Edney (12) when either of the two spring 
types acted alone. However, when both translational and 
rotational boundary springs were active, the two solutions 
of this study gave significantly lower frequencies. For all 
subsequent comparisons to published results for cases when 
restrained transverse motion was allowed at the edges, lower 
frequencies were again predicted by the methods used herein. 
Assuming no numerical error was present, since double checks 
uncovered none, this indicates that the choice of displacement 
functions becomes increasingly important. By considering the 
boundary conditions along with the choice of functions used in 
reference (12), the following is noted. For cases where the 
boundary conditions are between simply-supported and clamped 
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(ie. top curve of figures (6) and (7)), Warburton and Edney 
utilized a weighted sum of simply-supported and clamped beam 
functions. As sucn, the zero displacement edge condition was 
satisfied exactly while the moment condition is approximated 
through the range of elastic values with good results. On the 
lowest curve of those figures, a weignted sum of free and 
simply-supported beam functions were used. These functions 
satisfy the zero moment boundary condition exactly and 
approximate the shear condition over the range of elastic 
values. Again, by satisfying at least one of the boundary 
conditions, good results were achieved. For the intermediate 
curve where both types of spring were simultaneously active, 
the functions used in reference (12) were a weighted sum of 
simply-supported, clamped, free and sliding beam functions 
and collectively satisfy neither the shear or moment edge 
conditions. Rayleigh's principle guarantees the frequency 
predictions will be high by the argument of upper bounded 
eigenvalues. The weighted approximation on both types of 
elastic edge condition gives inflated results to the point 
of unpracticle accuracy. Unfortunately, one of the purposes 
of the published article was to give quick and easy frequency 
predictions with suitable accuracy. 

The additional comparisons given by tables (3), (4), (6) 
and (7) for cases when restrained transverse motion is allowed 
at the edges showed similar results. The solutions of this 
study gave lower predictions of frequency for all cases 
except those of table (6) where the finite element predictions 
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were inconsistent. The reason for the inflated results was 
discovered. Accidentally, the pair of parallel edges that 
varied from simply-supported to free states had been over 
constrained. For plates with unequal translational restraint 
on boundaries, transverse edge displacement is not equal on 
all edges. As a result, the boundary nodal displacements on 
the finite element model must have freedom of rotation about 

J 

an axis perpendicular to tne edge. This rotational freedom 
was inadvertently constrained for the computer runs of table 
(6). When this was corrected, the NASTRAN solution agreed 
with the elastic Ritz predictions. At a value of K*=10 in 
the top row of the table, the finite element dimensionless 
frequency parameter was found to be 11.25 compared to 11.31 
for the elastic Ritz solution. The SS-Fr-SS-Fr case has an 
exact value as given by Leissa. The finite element prediction 
with corrections returned a value of 9.659 as compared to 
9.676 and 9.6314 for the elastic Ritz and exact values, 
respectively. Therefore, the earlier major discrepancy of 
frequency prediction between NASTRAN and the elastic Ritz 
solutions may be attributed to user error. The minor dif- 
ferences of results between these two solutions can be 
attributed to convergence considerations. It is noted that 
the frequency results from the two approximations used in 
this study are not consistently nigher or lower in comparison. 
The single term Ritz solution is also guaranteed to be high 
by Rayleigh's principle while an insufficient number of 
elements on the finite element model will slow it's 
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convergence. For any given comparison, it is difficult to 
predict which of these factors will take precedence to give 
more agreeable estimates of frequency. 

Extension of the single term Ritz solution to predict 
response is generally less reliable than frequency prediction. 
Warburton (23) demonstrated that the single term solution may 
give very good results for displacement based on the conver- 
gence rate. The stress/strain convergence is usually slower 
however, especially for plates exhibiting unsvmmetric boundary 
conditions or aspect ratios not equal to one. In general, if 
convergence of natural frequencies is slow, the convergence 
for response calculations will be slower still. Therefore, 
for cases deviating far from the square plate with equal or 
symmetric edge conditions, a multiple mode calculation with 
more terms in the Ritz series is recommended. 

The fundamental resonance displacement and acceleration 
results showed good comparison between the single mode Ritz 
and finite element solutions in figures (8), (9), (10) and 
(11). This indicates that the convergence for displacements 
in the Ritz solution is rapid enough for good estimates. 
Combining good natural frequency and displacement comparisons, 
the accelerations predicted by both methods are expected to 
agree. Figures (10) and (11) show general agreement. 

The strain response values as seen in figures (12) and 
(13) are not as agreeable for the two methods in all cases. 
Several factors may account for this. The convergence of the 
single term Ritz solution is not guaranteed, as mentioned. As 
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well, the finite element solution results are only as good as 
the ref inementof the plate modeling and the element disp- 
lacement functions chosen. The deviation of NASTRAN strain 
values from the Ritz solution is seen to occur as restraining 
values fall below 10 3 on figures (12B) and (120. In these 
regions, the displacements begin increasing at a greater rate 
(see figure (8)). The continuous distribution of strain over 
the plate domain can only be approximated by the finite 
element model and since stress is constant over a given 
element, the distribution from NASTRAN is piecewise linear. 

As a result, for larger displacements, this piecewise dist- 
ribution will show greater deviation and is less reliable to 
approximate the continuous distribution. By increasing the 
number of elements or using mesh refinement in areas of 
interest, a better approximation should result. 

Another factor to help explain the disagreement of 
strain between the two solutions concerns the modeling of 
edge bending. For any given element, the corners are subject 
to differing magnitudes of translational and rotational disp- 
lacement. These differences result in shear and moment forces 
at the nodes. In the interior of the plate, continuity of 
displacement at the nodes joining elements cancels out these 
forces. That is, the forces on one element are equal and 
opposite to those of an adjoining element. At the edges, these 
effects are elastically restrained. As the stiffness of the 
elastic restraints are lessened, the forces at the edges of 
the model are neglected. As a result the edge elements are 
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subject to little or no bending. To properly account for this, 
it is suggested that numerically appropriate artificial support 
reactions and bending moments be applied at the boundaries to 
induce bending in the edge elements. This would increase the 
strain values returned by NASTRAN. 

Table (12) gave displacement and strain values for 
various combinations of edge restraint. The agreement was 
good. The NASTRAN finite element results were taken from 
reference (5) and non-dimensionalized. In that article, the 
edge strains were shown to increase as the transverse boundary 
spring stiffness was decreased and the rotational edge flexi- 
bility was increased. Since the experimental strains were 
lower than those predicted by a factor of three, the introduc- 
tion of restrained edge translation does not account for the 
bias error. The elastic Ritz results of this study have subs- 
tantiated that claim. 

For practical applications of the two approaches of 
this study, the following factors to be considered are sum- 
marized. The Ritz solution should be extended to a multiple 
mode solution with a qreater number of terms. Although, the 
frequencies from a single term solution can be accurate, tng 
convergence of solution to predict response will be enhanced. 
This will add to the complexity of solution since the off- 
diagonal terms in the frequency determinant and contributions 
to response of higher modes in the determination of amplitude 
coefficients will be included in the analysis. The use of 
subroutines available in math and science software libraries 
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to solve the eigenvalue problem and matrix simultaneous 
equations should ease the burden. Also, the finite element 
model should be adjusted with a greater number of elements 
or mesh refinement in the areas of stress calculations. An 
increased number of elements would allow for a better 
approximation of the mode shapes by the finite element model. 
The edges should also be modeled with proper reactions. 
Comments on the Use of. MACSYMA and NAS TRAN 

MACSYMA's evolution has been steady and methodical. The 
capabilities are continuously being expanded so that the 
program size is increasing rapidly. There are many problems 
in mathematics and in physical sciences to which MACSYMA has 
not been applied. For these reasons, it is impossible for one 
person to fully understand all of MACSYMA's capabilities. To 
the novice user, MACSYMA is somewhat of a "black box" that 
miraculously returns answers to complicated problems in math- 
ematics. 

In this study, MACSYMA was used to solve the simul- 
taneous equations of boundary conditions in order to specify 
the functions that describe the normal modes of vibration of 
the elastically restrained beam. This is the most difficult 
part of the solution. It is assumed that this difficulty is 
the reason no published material on the use of these beam 
functions in plate vibrations was found. Although the deriv- 
ation of the functions is difficult, it is not impossible to 
derive these functions by hand. In fact, they were derived by 
hand in order to assess that part of the contribution from 
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MACSYMA . The solution, when done by hand, comprised 14 pages. 
Care was taken to avoid mistakes (sign errors are a specialty) 
so that with simplifications and verification, the solution to 
get the beam functions took approximately 25 hours. When 
MACSYMA was applied, the most time was spent on problem setup. 
Once the input of the boundary equations was done and the 
general expressions for tne beam functions were substituted 
into those equations, MACSYMA returned the- completed solution 
for coefficients and the transcendental equation to determine 
the argument, alph(m), in about five minutes. Indeed, the 
expressions given by MACSYMA before simplification were 
fairly compact and programmable. 

MACSYMA was also used to perform all integration 
necessary to form the Rayleigh quotient and the work of the 
pressure load on the plate. The expressions from equations 
(22) through (26) for displacement and strain were derived 
completely on the computer. The entire analytical solution 
was done using MACSYMA. 

Once expressions were simplified to the desired degree, 
the Fortran code was generated. Fortran routines were written 
and debugged as the solution derivation progressed. Wnen the 
solution for the coefficients and argument of the beam func- 
tions was complete, the Newton-Raphson subroutine was written 
and verified. The necessary Fortran subroutine for calculation 
of the integral expressions and the formation of the Rayleigh 
quotient was written once progress of- analytic derivation per- . 
mitted. Finally the expressions for displacement and strain 



Brewer 83 


were derived and the response routine was programmed. The 
solution was complete. The derivation and programming of the 
entire Ritz algorithm was done on the computer. The major role 
of this user was to keep track of the necessary steps and to 
debug the Fortran program. 

One of the features of MACSYMA is the share directory. 

The share directory is a library of special math routines that 
were written by MACSYMA users and incorporated into the soft- 
ware for general use. As more physicists and engineers become 
aware of symbolic manipulation systems, tnese share direc- 
tories will begin to branch out. General solution algorithms 
to problems in the physical sciences will become available. 

For example, the Ritz procedure for plate vibrations might be 
incorporated into the system. The interactive user will input 
the functions and specify the boundary conditions. The soft- 
ware will grind out the necessary expressions and return a 
complete solution for natural frequencies. At the touch of a 
button, the Fortran code is given and numerical evaluations 
begin. Of course this is an over simplification but it is not 
an unrealistic concept. 

NASTRAN, like MACSYMA, is a large software package that 
no one user can completely grasp. It is therefore also a "black 
box" system. The user must be sure to follow the rules and 
format the input correctly in order to properly model the prob- 
lem. The results of table (6) were not corrected in this manu- 
script since they serve as an example. In the long formatted 
input listing, a single digit caused erroneous results. 
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One of the problems encountered in the course of this 
study was the length of time required and amount of storage 
needed by NASTRAN . Typical runs of NASTRAN, whicn required the 
loading of tne 15 modules that comprise the extensive package, 
used about five minutes of CPU time on the Cyber 175 system. 

In one instance, the number of degrees of freedom on the plate 
were increased nine-fold. That is, the number of elements was 
increased to check frequency convergence. The required CPU 
time was likewise increased nine-fold. Since the NASTRAN runs 
were made on a time share computer system, generating results 
was often frustrating and always time consuming. Runs were 
required to determine natural frequencies and again to apply 
an harmonic load at the fundamental frequency. As a result, 
NASTRAN did not lend itself conveniently to the parametric 
and response analyses of this study. 

An advantage of NASTRAN is that it can analyze complex 
problems. The finite element method provides a means to examine 
problems where no analytical method is available. The capabil- 
ities of NASTRAN are extensive but they have a limit. Finite 
element formulations to problems that NASTRAN can not presently 
handle are being derived by researchers everyday. It is only a 
matter of time that these solutions will be incorporated into 
software packages such as NASTRAN. 

A small fraction of the capabilities of both MACSYMA 
and NASTRAN were required for the analysis of this study. 
Therefore, this user must be considered a novice on both 
systems. The work did serve as a good introduction to the 
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available software. This introduction instilled an appreciation 
and respect for the capabilities of MACSYMA and NASTRAN that 
allows the following conclusions. 

Any engineer that is involved in complicated derivations 
of solution should become aware of symbolic manipulation 
systems and consider their use. The time saving alone is a 
valuable asset. These systems also eliminate the chances of 
mistake that accompany long algebraic and calculus related 
manipulations. From the point of view of this user, computer 
aided mathematics is tremendously powerful tool for the 
engineer . 

Finite element software has come of age. The capab- 
ilities of the available packages are always increasing. Even 
for problems where analytical solutions exist and are used, 
a finite element solution is a valid approach to check results 
or compare with experimental data. 
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APPENDIX A 
Notation 


It is notea here that MACSYMA does not recognize Greek 
nomenclature. Therefore, english abbreviations were used for 
those characters in all MACSYMA displays contained in this 
manuscript. Also, MACSYMA does not allow the use of upper 
and lower case characters simultaneously so that all MACSYMA 


displays are in lower case. For these reasons, the notation 
list below is expanded to include the MACSYMA notation used. 
The MACSYMA notations are enclosed in brackets '[]'. 


x,y, z 
a , b 

h 

t 


U), f 

• nf nuJ 
p , l rho) 


-z 


pL r ( pr ] 
n>, [d] 

w( x, y, t ) 


Cooridinate directions. 

Length (x direction) and width (y direction) 
of the plate, respectively. 

Thickness (z direction) of the plate. 

Time. 

Radian and cyclic frequencies, respectively. 

Density (mass/unit volume). 

Damping coeficient and % critical damping 
coeficient, respectively. 


Poisson's ratio. 


Flexural stiffness constant, D=Eh 3 /12(l -pL*) t 
where E==Young's Modulus. 

Transverse displacement (z direction) where 
W(x,y) is seperated spatial variation. 


-u.2 
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p ( x , y , t ) 


V,T 


Q 


Normally incident pressure on plate surface 
where P(x,y) = P is uniform spatial variation 
over the surface. 

Potential or strain energy and kinetic energy 
of tne plate, respectively. V is associated 
with elastic energy in the plate and boundary 
springs. 

Work done on the plate by pressure, p(x,y,t). 


X m (x),Y n (y) Beam functions in respective x and y directions 
[x(m), xm(x)j chosen for transverse displacement variation, 

( y (n ) , yn(y ) ) W(x,y ) . 


A mn ( a ( m , n ) ) 

B„, , C m , D_ Amplitude coeficients found in beam 

[S{mT!cTm) ,d(m) 1 

E n ,F n ,G n functions representing W(x,y). 

( e ( n) , f ( n ) ,g (n ) J 


am' Ulph(m) ] 
P„,l bet(n)) 

Ri , Ki 
[ri,ki ] 


Characteristic values found in the arguments 
of beam functions. 

Rotational and transverse translational 
restraint values, respectively. 


r,i,m, Indices, (positive integers). 

n,p,q 


<j , g Stress and strain, respectively. 

( cc,ss ,chch 

shsh,cs,cch Integrals defined in Appendix C. 

csh,sch,ssn 

chsh | 


[ xhm,x jm,xkm 

xlm, xosq , xasq Integrals of beam functions defined in 
xposq,xpasq Appendix C. 

xqm | 

[%e, %pil Euler constant, e-2 . 17182818 and 

TT =3. 14159265, respectively. Recognized 
by MACSYMA as given. 



Biharmonic differential operator defined as. 



0X 4 (3 x 2 y 2 (3 y 4 
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APPENDIX B 

Derivation of Beam Functions 


The equations defining the elastic boundary conditions 
that the assumed displcicement function, W(x,y) must satisfy 
were given in equations (6', 7'). These equations may be 
seperated so that X m (x} must satisfy equations (6A', 6B' 

7A ' , 7B') and Y n (y) must satisfy equations (6C' , 6D', 7C ' , 
7D'). Once the coefficients, (B, C, D) and argument, alph(m), 
of X m are obtained, the corresponding coefficients, (E, F, 

G) and argument, bet(n), of Y n can be determined by the 
substitution described in chapter 2. 

Rewriting the boundary condition equations for Xjjj and 
utilizing MACSYMA we get, 

2 

d d 

( — (xm(o) ) ) rl * d ( (xm(o))) ( 31 ) 

dx 2 

dx 

2 

d d 

(-- (xzn(a) ) ) r2 « - { ( xm ( a ) ) ) d 

dx 2 (B2) 

dx 

3 
d 

kl xm(o) * - d ( (xm(o) ) ) 

3 

dx 


(B3) 
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xro(a) k2 


3 



3 


dx 


(xm(a) ) ) d 


(B4 ) 


Substituting the expression for given in equation ( 10 ), 
differentiating and evaluating the appropriate equations at 
x=o and x=a, equations (Bl, B2, B3, B4) become. 


alph(m) c(m) alph(m) b(m) alph (m) d(m) 

(~ + ■—) rl = d ( — 

a a 2 

a 


alph (m) 

2 

a 

( B5 ) 


alph(m) d (m) sinh(alph(m) ) alpn(m) sin(alph(m)) 

(-- * - 

a a 

alph(m) b(m) cosh(alph(m) ) alph(m) c(m) cos(alph(m) ) 

— . + — 

a a 


) r2 


- d ( 


2 2 
alph (m) b(m) sinh(alph(m) ) alph (m) c(m) sin(alph(m)) 


2 2 
alph (m) d(m) cosh (alph(m) ) alph (m) cos (alph(m) ) ^ 

) 

2 2 
a a 


(B6) 


3 3 

alph (m) b(m) alph (m) c(m) 

kl (d(m) + 1) = - d ( ) 

3 3 
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k2 (b(m) sinh(alph(m) ) + c(m) sin(alph(m)) 
+ d(m) cosh(alph(m) ) ♦ cos (alph(m) ) ) » 


3 

alph (m) d(m) sinh( alph(m) ) 
d ( 


3 


a 


3 

alph (m) sin(alph(m)) 
+ — - — 

3 

a 


3 

alph (m) b(m) cosh(alph(m) ) 
+ — -- — ........ — -- 

3 

a 


3 

alph (m) c(m) cos(alph(m)) 


3 


) 


a 


(B8 ) 

Using equations (B5, B7, B8) and solving simultaneously for 
B m , C m and D m , MACSYMA gives. 


4 3 

b(m) * (- (a d kl alph (m) 

(sinh(alph(m) ) - sin(alph(m) ) ) 

7 

+ a kl (k2 cos(alphtm)) - k2 cosh( alph(m) ) ) ) rl 
3 7 

- d alph (m) (sinh( alph(m) ) + sin( alph(m) ) ) 

6 

+ 2 a d kl k2 alph(m) sin(alph(m)) 

3 2 4 

-a a alph (m) {- k2 cosh (alph (m) ) - k2 cos(alph(m)) 

- 2 kl cos(alph(m) ) ) ) /denominator 

(B9) 
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4 3 

c(m) = ( (a d kl alph (m) 

(sinh(alph(m) ) - sin(alph(m) ) ) 

7 

+ a kl (k2 cos(alph(m) ) - k2 cosh ( alph (m) )) ) rl 
6 

-2a d kl k2 alph(m) sinh(alph(m) ) 

3 7 

+ d alph (m) (- sinh(alph(m) ) - sin(alph(m) ) ) 

3 2 4 

♦ a d alph (m) (k2 cos*n(alph(m) ) + 2 kl cosh( alph(m) ) 

+ k2 cos (alpn(m) ) )) /denominator (BIO) 


7 

d (m) * (- (a kl (k2 sinh(alph(m) ) - k2 sin( alph(m) ) ) 
2 6 

+ 2 a d alph (m) sin(alph(m) ) 

4 3 

+ a d alph (m) (kl ( - cosh(alph(m) ) - cos(alph(m) ) ) 

3 2 4 

- 2 k2 cost alph (m) ) ) ) rl - a d alph (m) 

(k2 sinh(alphtm) ) + k2 sin(alph(m) ) ) 

3 7 

- d alph (m) (cos(alph(m) ) - cosh(alph(ro) ) ) ) /denominator 


(Bll) 

where, 


7 

denominator * (a kl (k2 sxnh(alph(m) ) 

2 6 

- k2 sin ( alph (m) ) ) ♦ 2 a d alph (m) sinh(alph(m) ) 


+ a d alph (m) (kl (- cosh(alph(m) ) - cos(alph(m) ) ) 
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3 2 4 

- 2 K2 cosh(alph(rn) ) ) ) rl + a d alph (in) 

( - k2 sinh(alph(m> ) - k2 sin(alph(m) ) ) 

3 7 

♦ d alph (m) (cosh(alph(m) ) - cos ( alph (m) ) ) (B12) 

Equation (B6) can be rewritten as, 

alph(m) d(m) sinn(alph(m) ) alph(m) sin( alph(m) ) 

( 

a a 

alph(m) b(m) cosh(alph(m) ) alph(m) c(m) cos(alpn(m) ) 

+ — + •) 

a a 

2 

alph (m) b(m) sinh(alph(m) ) 

r2 + d ( 

2 

a 

2 2 
alph (m) c(m) sin(alph(m) ) alph (m) d(m) cosh(alph(m) ) 

— — — + — — 

2 2 

a a 

2 

alph (m) cost alph ( in ) ) 

) « 0 

2 

a 

(B13) 

and substitution of equations (B9, BlO, Bll, B12) into (B13) 
results in the transcendental equation that is used to solve 
for alph(m). 

The solution of equation (B13) was done numerically by 
Newton-Raphson method which necessitates the determination 
of the derivative of the left-hand side of equation (B13) 
with respect to alph(m). MACSYMA was used for this calculation, 
however the expression is quite long and not included here. 
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APPENDIX C 

Calculation of Beam Integrals 

The integrals necessary in the calculation of the 
Rayleigh quotient for natural frequencies and subsequently 
the work done on the plate for forced response were given 
in equations (2, 3, 16, 17 and 19). These integrals are 
evaluated by inserting the expressions assumed for the mode 
displacement functions (equations (10 and 11). Since these 
functions are separable and similar in x and y, the integrals 
may be evaluated independently for x and similar expressions 
for y are obtained by the substitution described in chapter 2. 

Inserting the expression for (x) into the integrands 
of the appropriate integrals and expanding these integrand 
expressions, the following integrals are needed. 



(C2) 
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a 

/ 

[ 2 alph(m) x 

chch » I cosh ( ) dx 

1 a 

/ 

0 (C3) 


/ 

{ 2 alph(ro) x 

shsn * I sinh ( ) dx 

] a 

/ 


( C4 ) 


/ 

( alph(m) x alph(m) x 

cs = I cos ( ) sin( ) dx 

J a a 

/ 

0 


(C5) 


/ 

( alph(m) x alph(m) x 

ccn = I cos ( ) cosh ( ) dx 

] a a 

/ 

0 


<C6) 


/ 

[ alph(m) x alph(m) x 

csh * I cost ) sinh ( ) dx 

] a a 

/ 

0 


( C7 ) 


/ 

[ alph(m) x alph(m) x 

sch * I cosh ( ) sin ( ) dx 

1 a a 

/ 

0 


(C8 ) 
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ssh = 


/ 

[ 

I 

J 

/ 


alph(m) x alph(m) x 

s i n ( ~ - — ■ — — ) sinh ( ) dx 

a a 


( C9 ) 


/ 

l alph(m) x alph(m) x 

chsh = I cosh ( ) sinh ( ) dx 

] a a 

/ 

0 (CIO) 

In general, MACSYMA converts the hyperbolic sine and 

hyperbolic cosine (sinh and cosh) into exponential form as 

an aid to integration by the familiar identities, 

U - u u — u 

%e + %e %e - %e 

cosh(u) * rsinh(u) 

2 2 ( Cl 1 ,C1 2 ) 

After integration, MACSYMA does not recombine the exponential 
terms and gives. 


cc 


a sin( 2 alph(m)) + 2 a alpn(m) 
4 alph(m) 


(Cl 3 ) 


ss = 


a sin( 2 alph(m) ) - 2 a alph(m) 
4 alph(m) 


(C14 ) 


te 


- 2 alph(m) 4 alph(m) 2 alph(m) 

(a %e + 4 a alph(m) %e - a) 


chch 


8 alph(m) 


(Cl 5) 



Brewer 99 


shsh 


- 2 alph(m) 4 alph(m) 2 alpMm) 

%e (a %e - 4 a alph(m) %c 


8 alph(m) 


- a 


( Cl 6 ) 


cs * 


a sin (alph(m) ) 
2 alph(m) 


(Cl 7) 


alph(m) - alph(m) 

a %e sin( alph(m) ) a %e sin(alph(m)) 

cch = + 

4 alph(m) 4 alph(m) 

alph(m) -alph(m) 

a %e cos (alph(m) ) a %e cos(alph(m) ) 


4 alph(m) 


4 alph(m) 


(C18) 


csh = 


alph(m) - alph(m) 

a %e sin( alph(m) ) a %e sin(alpn(m)) 


a %e 


4 alph(m) 
alph(m) 


4 alph(m) 
- alph(m) 


cos(alph(m) ) a %e 



4 alph(m) 4 alph(m) 


cos (alph(m) ) 


2 alph(m) 
alph(m) 


a %e 


sch = 


(Cl 9) 

- alph(m) 

sin(alph(m) ) a %e .sin(alph{m) ) 


a %e 


4 alph(m) 
alph(m) 


4 alph(m) 


- alph(m) 

cosCalph(m) ) a %e cos (alpn(m) ) 


4 alph(m) 


4 alph(m) 


a > 


2 alph(m) 


(C20) 
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alph(m) - alph(m) 

a %e sin(alph(m)) a %e sin(alph(m) ) 

ssh — + * 

4 alph(m) 4 alph(m) 

alph(m) - alph(m) 

a %e cos (alph(m) ) a %e cos(alph(m) ) 

- — — — — — — 

4 alph(m) 4 alph(m) 

(C21 ) 

2 

a cosh (alph(m) ) a 

chsn ■ 

2 alph(m) 2 alph(m) 

(C22 ) 


By recognizing the hyperbolic terms and simplifying 
according to the identities of equations (Cll and Cl 2) the 


exponential expressions may be reduced to, 

chch = (a/4alph(m)) ( sinh( 2alph(m) ) +2alph(m) ) (C23) 

shsh = (a/4alph(m) ) ( sinh( 2alph(m) )-2alph(m) ) (C24) 

cch = (a/2alph(m) ) (sin(alph(m) )cosh(alph(m) ) 

+cos(alph(m) )sinh(alph(m) ) ) (C25) 

csh = (a/2alph(m)) (sin(alph(m) )sinh(alph(m) ) 

+cos(alpn(m) )cosh(alph(m) ) — 1 ) (C26) 

sch = (a/2alph(m)) (sin(alph(m) )sinh( alph(m) ) 

-cos(alph(m) )cosh(alph(m) )+l ) (C27) 

ssh = (a/2alph(m)) (sin(alph(m) )cosh(alph(m) } 

-cos(alph(m) ) sinh(alph(m) ) ) (C28) 


The integrals that result from the insertion of X m (x) 
are then defined as. 
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/ 

l 2 

xhm * I xm (x) dx 

] 

/ 

0 


/ 

C d 2 

xjm * I ( — (xm(x) ) ) dx 
I dx 

/ 

0 


a 

/ 2 

l d 2 

xkm = I ( ( xm ( x ) ) ) dx 

) 2 
/ dx 
0 


/ 2 

[ d 

xlm * I xm(x) ( (xm(x))) dx 

) 2 
/ dx 

0 


and can be shown to be given by, 
xhm = b(m) b(m) shsh + 2 b(m) c(m) ssh + 2 b(m) d(m) 
+ 2 b(m) csh + c(m) c(m) ss + 2 c(m) d(m) sch 
+2 c(m) cs + d(m) d(m) chch + 2 d(m) cch + cc 
2 

xjm = (alph(m)/a) (d(m) d(m) shsh - 2 d(m) ssh 
+ 2 b(m) d(m) shch + 2 c(m) d(m) csh + ss 
- 2 b(m) sch - 2 c(m) cs + b(m) b(m) chch 
+ 2 b(m) c(m) cch + c(m) c(m) cc) 


(C29 ) 
(C30) 

(C31 ) 

(C32) 

shch 

(C33) 


(C34) 



Brewer 1 02 


4 

xkm = (alph(m)/a) (b(m) b(m) shsh - 2 b(m) c(m) ssh 

+ 2 b(m) d(m) shch - 2 b(m) csh + c(m) c(m) ss 

- 2 c(m) d(m) sch + 2 c(m) cs + d(m) d(m) chch 

- 2 d(m) cch + cc) 

2 ( C35 ) 

xlm = (alph(m)/a ) (b(m) b(m) shsh + 2 b(m) d(m) shch 

- c(m) c(m) ss - 2 c(m) cs + d(m) d{m) chch - cc 

(C36 ) 

The final necessary expressions for the determination of 
natural frequencies arise from the elastic restraint energy 
terms of equations (16 and 17). Letting [X m (o) J^xosq, 

(X m (a) ] 2 =xasq, [ X ' m ( o) ] 2 =xposq and [ X ' m ( a ) ] 2 =xpasq where 
indicates the derivative with respect to x, then. 


2 

xosq * (d(m) ♦ 1) 

2 2 

xasq = b (m) sinh (alph(m) ) 

♦ 2 b(m) c(m) sin(alph(m) ) sinh(alpn(m) ) 

+ 2 b(m) d(m) cosh(alph(m) ) sinh(alph(m) ) 


(C37 ) 


2 2 

+ 2 b(m ) cos (alph(m) ) sinh(alph(m) ) + c (m) sin (alph(m)) 

+ 2 c(m) d(m) cosh( alph(m) ) sin(alph(m) ) 

2 2 

+ 2 c(m) cos(alph(m) ) sin(alpn(m) ) + d (in) cosh (alph(m)) 

2 

♦ 2 d ( m ) cos( alph(m) ) cosh(alph(m) ) + cos (alph(m) ) (C38) 


2 2 
alph (m) (c(m) + b(ni)) 


xposq 


2 
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2 2 2 

xpasq = alph (m) (d (m) sinh (alph(m)) 

* 2 d(m) sin(alph(m) ) sinh(alph(m) ) 

+ 2 b(m) d(m) cosh(alph(m) ) sinh(alph(m) ) 

2 

+ 2 c(m) d(m) cos (alph(m) ) sinh(alpn(m) ) + sin (alph(m) ) 

- 2 b(m) cosh(alph(m) ) sin(alph(m) ) 

2 2 

- 2 c(n>) cos (alph(m) ) sin(alph(m) ) + b (m) cosh (alph(m)) 

+ 2 b(m) c(m) cos (alph(m) ) cosh(alph(m) ) 

2 2 2 
+ c (m) cos (alph(m]i ) )/a 

(C40) 


Again, similar expressions may be obtained for integrals 
of Y n (y ) and it's derivatives by the substitution described 
on page 18 (chapter 2]i. 

The Rayleigh quotient can be formed as, 

2 

(j = (d (xkm yhn + xhm ykn + 2 f± xlm yin + 2 (1-^i,) xjra yjn) 

+ Rj xposq yhn + R 2 xpasq yhn + R 3 yposq xhm + R 4 ypbsq xhm 

+ Kj xosq yhn + K 2 xasq yhn + K 3 yosq xhm + K 4 ybsq xhm ] ^ 

Ip h xhm yhn ] 

r (C41) 


Finally, the integration defined by equation (19) for 
the work of the uniform pressure on the plate is needed in 
order to calculate the amplitude coeficients, A mn . This 
integral can be written for X m (x) as, 

a 

/ 

( 

xqm » I xm( x) dx 
1 
/ 

0 


(C42) 
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and is evaluated as, 

xqm * a (d(m) sinh( alph(m) ) ♦ sin(alpMm) ) 

♦ b(m) cosh(alph(m) ) - c(m) cos (alph(m) ) + c(m) - b(m)) 
/alph(m) 


(C43) 
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APPENDIX D 

Solution for Coeficients of Rayleigh-Ritz Series 


The equation used to determine the coeficients (A^) of 
the Rayleigh-Ritz series was given iq chapter 2 (equation 
(22)). Substitution of the assumed displacement functions 
(equations (10 and 11) and the integral expressions from 
Appendix C, equation (22) can be written as. 





Ajnn xqm vqn) 

(Dl) 


or , 


2 2 

(Ajjjjj p hW xhm yhn)(l - Jj^nn) = -P xqm yqn 


Solving for A jnn , 
\in ~ 


00 


-P xqm yqn 


(D2) 


2 /. 2 

PhCO xhm yhn (1 - %in ) 


U) 


-P xqm yqn^ 
P h xhm yhn ( ^mn 


T 
W ) 


(D3) 


Equation (D3) is for the undamped case. Damping can be 
incorporated into the solution in the form of percent critical 
damping by adding i2/*CJ mn CU to tWmn -V ) in the denominator. 


